maptrian.cpp

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

//************************************************************************
//maptrian - partial update of lower triangular factor for a sym. matrix
//************************************************************************
#include <malib.h>

//************************************************************************
int EmsMatrix::maptrian(void)
{
   int rc = SUCCESS;     
   int npartial = 0;                 
   int i;                      
   struct mastatest *mastatept;     
   struct maorderst *maorderpt;      
  
   for (i = 1; i<= nmastate;i++)
   {
      mastatept=mastate_vt[i-1];
	  if ( mastatept->stat.partial )
      {
         mastatept->stat.partial = FALSE;
         maorderpt = maorder_vt[mastatept->maorder_itm - 1];
         if ( maorderpt->stat.active )
   	     maorderpt->stat.partial = TRUE;
      }
   }
   
  
   for (i=0;i<nmastate;i++)
   {

	  maorderpt = maorder_vt[i];
      if ( maorderpt->stat.partial )
      {
         npartial++;
         rc = lfinitialize_row (maorderpt);
         if ( rc != SUCCESS )
   	 break;
   
         rc = lfevaluate_row (maorderpt);
         if ( rc != SUCCESS )
   	 break;
   
         rc = lfnormalize_row ( maorderpt);
         if ( rc != SUCCESS )
   	 {
            printf("maptrian: trianggulizing G-matrix encountered single matrix\n");
   	    break;
   	 }
      }
   }
   
   return ( rc );
}

//************************************************************************
// initialize_row - initialize working row from C-matrix coefficients
//************************************************************************
int EmsMatrix::lfinitialize_row(struct maorderst *tmaorderpt)
{
   int rc = SUCCESS;           
   int i, j;                      
   int cmtrxcnt, magmtrx_cnt;               
   bool error;              
   struct maamtrxst *cmtrxpt; 
   struct magmtrxst *magmtrxpt, *magmtrx_pt;
   struct maorderst *maorder_pt; 
   
   cmtrxcnt = tmaorderpt->cmtrxcnt;
   magmtrx_cnt = tmaorderpt->magmtrx_cnt;
  
   for (j = 0;j<magmtrx_cnt;j++)
   {    
       magmtrxpt = magmtrx_vt[j+tmaorderpt->magmtrx_itm - 1];	
	   magmtrxpt->element = 0.0;
   }

   cmtrxpt = (struct maamtrxst *)cmtrxft + tmaorderpt->cmtrxitm - 1;

   for (i = 0, j = 0; i < cmtrxcnt; i++, cmtrxpt++ )
   {
 	  maorder_pt = maorder_vt[cmtrxpt->colitm-1];
      if ( maorder_pt->stat.active )
      {
     
	  for ( error = TRUE; j < magmtrx_cnt; j++ )//magmtrxpt++
   	  {
            magmtrxpt = magmtrx_vt[j+tmaorderpt->magmtrx_itm - 1];;
		    
			if ( cmtrxpt->colitm == magmtrxpt->colitm )
            {
               magmtrxpt->element = cmtrxpt->element;
               error = FALSE;
               j++;
               break;
            }
   	 }
   
         if ( error )
   	 {
            rc = AS_TRI_CMATRIX_CHG;
            printf("maptrian: tri. Cmatrix - Cmatrix struc. changed \n" );
            break;
   	 }
      }
   }
   
   return ( rc );
}

//************************************************************************
// evaluate_row - evaluate working row elements
//************************************************************************
int EmsMatrix::lfevaluate_row (struct maorderst *tmaorderpt)
{
   int i, j;                      
   int rc = SUCCESS, row;                    
   int magmtrx_itm, magmtrx_cnt;               
   bool error = FALSE;       
   struct magmtrxst *magmtrxpt, *magmtrx_pt, *gpt;      
   struct maorderst *maorder_pt;  
   double product;             
   
   magmtrx_cnt = tmaorderpt->magmtrx_cnt;
   row = mastate_vt[tmaorderpt->mastate_itm - 1]->maorder_itm;
  
   for ( i = 0; i < magmtrx_cnt; i++)
   {
      magmtrxpt = magmtrx_vt[i + tmaorderpt->magmtrx_itm - 1];
	  if ( magmtrxpt->element != 0.0 )
      {
         gpt = magmtrxpt + 1;
         j = i + 1;
   
         magmtrx_itm = maorder_vt[magmtrxpt->colitm - 1]->magmtrx_col;
         while ( magmtrx_itm != 0 )
         {
            magmtrx_pt = magmtrx_vt[magmtrx_itm - 1];
            if ( magmtrx_pt->rowitm >= row )
               break;
           
            product = magmtrxpt->element * magmtrx_pt->element;
            for ( error = TRUE; j < magmtrx_cnt; j++, gpt++ )
            {
               if ( gpt->colitm == magmtrx_pt->rowitm )
               {
                  gpt->element -= product;
                  error = FALSE;
                  break;
               }
            }
   
            if ( error )
               break;
            magmtrx_itm = magmtrx_pt->magmtrx_itm;
         }
   
         lftag_row ( magmtrxpt);
   
         if ( error )
         {
            rc = AS_TRI_CMATRIX_CHG;
            printf("maptrian: tri. Cmatrix changed - Cmatrix structure changed \n");
            break;
         }
      }
   }
   
   if ( tmaorderpt->magmtrx_col != 0 )
   {
      magmtrxpt =  magmtrx_vt[tmaorderpt->magmtrx_col - 1];
      maorder_pt = maorder_vt[magmtrxpt->rowitm - 1];
      if ( maorder_pt->stat.active )
         maorder_pt->stat.partial = TRUE;
      else
         lftag_row ( magmtrxpt);
   }
   
   return ( rc );
}
   
//************************************************************************
// tag_row - tag rows on factor path
//************************************************************************
void EmsMatrix::lftag_row(struct magmtrxst *magmtrxpt)
{
   struct maorderst *maorderpt;  

   while ( magmtrxpt->magmtrx_itm  != 0 )
   {
      magmtrxpt = magmtrx_vt[magmtrxpt->magmtrx_itm - 1];
      maorderpt  = maorder_vt[magmtrxpt->rowitm - 1];
      if ( maorderpt->stat.active )
      {
         maorderpt->stat.partial = TRUE;
         break;
      }
   }
   return;
}

//************************************************************************
// normalize_row - normalize working row elements
//************************************************************************
int EmsMatrix::lfnormalize_row(struct maorderst *tmaorderpt)
{
   int rc = SUCCESS; 
   int i, magmtrx_cnt;                
   double norm, dot_product = 0.0;       
   struct magmtrxst *magmtrxpt;      
   struct maorderst  *maorder_pt;      
   struct mastatest *mastatept;      
   
   magmtrx_cnt = tmaorderpt->magmtrx_cnt;
   
   //magmtrxpt = maparmpt->magmtrxpt + maorderpt->magmtrx_itm - 1;
   
   for (i = 0; i < magmtrx_cnt; i++)
   {
      magmtrxpt = magmtrx_vt[i+tmaorderpt->magmtrx_itm - 1];
	  maorder_pt =maorder_vt[magmtrxpt->colitm - 1];
      norm    = magmtrxpt->element * maorder_pt->diag;
      dot_product += norm * magmtrxpt->element;
      magmtrxpt->element = norm;
   }
    
   tmaorderpt->stat.partial = FALSE;
   mastatept = mastate_vt[tmaorderpt->mastate_itm - 1];
   tmaorderpt->diag = mastatept->diag - dot_product;
   if ( fabs ( tmaorderpt->diag ) < AS_TRI_TOL )
   {
      rc = AS_TRI_SINGULAR;
      divg_state_pt = tmaorderpt->mastate_itm;
   }
   else
      tmaorderpt->diag = 1.0 / tmaorderpt->diag;
   
   return ( rc );
}