mainvert.cpp

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

//************************************************************************
// mainvert - sparse matrix inversion for non-zero elements in sym. matrix
// 矩阵转置
//************************************************************************
#include <malib.h>

//************************************************************************
int EmsMatrix::mainvert( void)
{
   int rc;    

   rc = check_parms();

   if ( rc == SUCCESS )
   {
      calc_matmtrx();

      proc_ainv_calc();
       
      if ( matmtrxft != NULL )
      {
	   free(matmtrxft);
	   matmtrxft = NULL;
      }
   }

   return(rc);
}

//************************************************************************
// check_parms - validate caller-initialized matrix parameters
//************************************************************************
int EmsMatrix::check_parms(void)
{
   int rc = SUCCESS; 
    
   if ( magmtrxft == NULL )
   {
      printf("mainvert: ERROR - magmtrxpt is NULL \n");
      rc = !SUCCESS;
   }

   if ( maorderft == NULL )
   {
      printf("mainvert: ERROR - maordrpt is NULL \n");
      rc = !SUCCESS;
   }

   if ( nmagmtrx == 0 )
   {
      printf("mainvert: ERROR - nmagmtrx is 0\n");
      rc = !SUCCESS;
   }
    
   if ( nmastate == 0)
   {
      printf("mainvert: ERROR - maparm.nmastatept is NULL \n");
      rc = !SUCCESS;
   }
  
   return(rc);
}

//************************************************************************
// calc_matmtrx - calculate elements of the T-matrix
//************************************************************************
void EmsMatrix::calc_matmtrx(void)
{
   int i,m_recNum;                         
   struct magmtrxst *magmtrxpt;     
   struct matmtrxst *matmtrxpt;     

   m_recNum = nmagmtrx;
   matmtrxft =(char *) new struct matmtrxst[m_recNum ];
   memset( matmtrxft, 0, m_recNum*sizeof(struct matmtrxst));
   matmtrxpt =(struct matmtrxst *)matmtrxft;
   for ( i = 0;i <magmtrx_vt.size( ); i++,matmtrxpt++ )
   {
      magmtrxpt = magmtrx_vt[i]; 
	  matmtrxpt->element  = -magmtrxpt->element;
      matmtrxpt->matmtrx_itm =  magmtrxpt->magmtrx_itm;
      matmtrxpt->row = magmtrxpt->colitm;
      matmtrxpt->col = magmtrxpt->rowitm;
   }
   return;
}

//************************************************************************
// proc_ainv_calc - process sparse A-inverse calculation
//************************************************************************
void EmsMatrix::proc_ainv_calc(void)
{
   register int i,j;        
   struct maorderst *tmaorderpt;
   double ainv_diag;      

   for (i =nmastate,j=0; i>0;i--,j++)
   {
	  tmaorderpt=maorder_vt[nmastate-1-j];
	  if ( tmaorderpt->stat.active )
      {
		calc_diag(tmaorderpt,&ainv_diag);
		proc_off_diags(i, tmaorderpt,ainv_diag);
      }
   }
   return;
}

//************************************************************************
// calc_diag - calculate A-inverse diagonal elements
//************************************************************************
void EmsMatrix::calc_diag(struct maorderst *tmaorderpt,double *ainv_diagpt)
{
   bool more_terms = TRUE; 
   struct magmtrxst *tmagmtrxpt;
   struct matmtrxst *tmatmtrxpt;
   struct matmtrxst *tbasept;     

   *ainv_diagpt = tmaorderpt->diag;
   tbasept = (struct matmtrxst *)matmtrxft;

   if ( tmaorderpt->magmtrx_col != 0 )
   {
      tmatmtrxpt = tbasept + tmaorderpt->magmtrx_col-1;
      tmagmtrxpt = magmtrx_vt[tmaorderpt->magmtrx_col - 1];

      while ( more_terms )
      {
	 if ( tmagmtrxpt->rowitm == tmatmtrxpt->col )
         {
	    *ainv_diagpt += tmatmtrxpt->element * tmagmtrxpt->element;
	    if ( tmatmtrxpt->matmtrx_itm != 0 )
	       tmatmtrxpt = tbasept + tmatmtrxpt->matmtrx_itm;
	    else
	       more_terms = FALSE;
         }
	 else if ( tmagmtrxpt->rowitm>tmatmtrxpt->col )
         {
	    if ( tmatmtrxpt->matmtrx_itm != 0 )
	       tmatmtrxpt = tbasept+tmatmtrxpt->matmtrx_itm-1;
	    else
	       more_terms = FALSE;
         }
	 else
         {
	    if ( tmagmtrxpt->magmtrx_itm != 0 )
	       tmagmtrxpt = magmtrx_vt[tmagmtrxpt->magmtrx_itm-1];
	    else
	       more_terms = FALSE;
         }
      }
   }

   tmaorderpt->diag = *ainv_diagpt;

   return;
}

//************************************************************************
// proc_off_diags - process calculation of A-inverse off-diagonal elements
//************************************************************************
void EmsMatrix::proc_off_diags(int current_i,struct maorderst *maorderpt,double ainv_diag)
{
   int i,j;                       
   struct magmtrxst *tmagmtrxpt;    

   for ( i = maorderpt->magmtrx_cnt,j=0 ; i > 0; i--,j++)
   { 
	   tmagmtrxpt = magmtrx_vt[maorderpt->magmtrx_itm + maorderpt->magmtrx_cnt-2-i];
	   calc_off_diag(current_i, maorderpt, tmagmtrxpt, ainv_diag);
   }
   return;
}

//************************************************************************
//calc_off_diag - calculate A-inverse off-diagonal element
//************************************************************************
void EmsMatrix::calc_off_diag(int current_j, struct maorderst *maorderpt, struct magmtrxst *zmtrxpt,double ainv_diag)
{
   int zi, zj;                       
   double off_diag = 0.0, z_element;             
   int magmtrx_itm;                 
   struct maorderst *wkmaorderpt;    
   bool more_terms = TRUE;   
   struct matmtrxst *matmtrxpt;
   struct matmtrxst *tbasept; 

   zi = zmtrxpt->colitm;
   zj = current_j;
   magmtrx_itm = maorderpt->magmtrx_itm;
   tbasept = (struct matmtrxst *)matmtrxft;

   if ( magmtrx_itm != 0 )
   {
      wkmaorderpt = maorder_vt[zi-1];
      matmtrxpt  = tbasept + wkmaorderpt->magmtrx_col-1;
       
      while ( more_terms )
      {
	 z_element = get_zkj(matmtrxpt->col, zj, magmtrx_itm,ainv_diag);
	 off_diag += matmtrxpt->element * z_element;
	 if ( matmtrxpt->matmtrx_itm != 0 )
	    matmtrxpt = tbasept + matmtrxpt->matmtrx_itm-1;
	 else
	    more_terms = FALSE;
      }
   }
   zmtrxpt->element = off_diag;

   return;
}

//************************************************************************
// get_zkj - get Z(k,j) term for current k and j
//************************************************************************
double EmsMatrix::get_zkj(int zk, int zj, int magmtrx_itm, double ainv_diag)
{
   struct maorderst *wkmaorderpt;     
   struct magmtrxst *wkmagmtrxpt;
   double z_element = 0.0;        
   int z,i;                

   if ( zk < zj )
   {
      wkmagmtrxpt = magmtrx_vt[magmtrx_itm - 1];
      z = zk;
      i=0;
      while ( wkmagmtrxpt->colitm < z )
	  {
		  i++;
		  //wkmagmtrxpt++;
		  wkmagmtrxpt=magmtrx_vt[magmtrx_itm - 1+i];
	  }
      if ( wkmagmtrxpt->colitm == z )
	 z_element = wkmagmtrxpt->element;
   }
   else if ( zk > zj )
   {
      wkmaorderpt = maorder_vt[zk - 1];
      if ( wkmaorderpt->magmtrx_itm != 0 )
      {
	 wkmagmtrxpt = magmtrx_vt[wkmaorderpt->magmtrx_itm - 1];
	 z = zj;
	 i=0;
	 while ( wkmagmtrxpt->colitm < z )
	 {
         i++;	
		 //wkmagmtrxpt++;
         wkmagmtrxpt = magmtrx_vt[i+wkmaorderpt->magmtrx_itm - 1];

	 }	  
	 if ( wkmagmtrxpt->colitm == z )
	    z_element = wkmagmtrxpt->element;
      }
   }
   else
      z_element = ainv_diag;

   return(z_element);
}