maoptim.cpp

来自「电力系统稀疏矩阵计算类,用法范例代码,包括三个目录,MALIB(类),MATES」· C++ 代码 · 共 796 行 · 第 1/2 页

//===========================================================================
//
//                            EMS高级应用软件
//
//===========================================================================
//  rev    0.0   7-29-2004      initial release   Qwbin
//  版本   0.0   7-29-2004      第一次发布        齐文斌
//         1.0   1-27-2006      修改为C++vector方式     齐文斌
//---------------------------------------------------------------------------

#include <malib.h>

//************************************************************************
// maoptim - sparse matrix optimal ordering
// 稀疏矩阵优化排序----主控 
// maorder_vt.push_back(maorderpt); 暂时放此
//************************************************************************
int EmsMatrix::maoptim()
{
   int rc;                       
   int n;                        
   int nonsymmetric;             
   struct opworkst work;     //临时空间,只有一个结构(一条纪录)     
   
   memset ( &work, 0, sizeof(work));

   opinitialize (&work);
   
   nonsymmetric = test_symmetry(&work);
   
   if ( nonsymmetric )
      make_symmetric ( &work );
   
   if ( ! work.stat.overflow )
   {
      scheme1 ( &work );
   
      while ( work.num_zero_fillin < work.num_zero_diagonal )//work.num_zero_fillin初始清零  mdn 8.11
      {
         n = work.num_zero_fillin;
         reorder_scheme1 ( &work );
   
         if ( work.num_zero_fillin == n )
         {
            printf("maoptim:REORDER FOR ZERO DIAGONAL FAILED %d \n", n );
            break;
         }
      }
      scheme2 ( &work );
   }

   opterminate ( &work );

   if ( work.stat.overflow )
      rc = !SUCCESS;
   else
      rc = SUCCESS;
   
   return ( rc );
}

//************************************************************************
// initialize - initialize internal work areas
// 初始化,给 workpt赋初值
// 统计每列连接元素数statept->connection_cnt
//************************************************************************
// 函数名称：       opinitialize 
// 返回值:          无
// 函数功能：       形成链表，建立statest与linkst的关联结构，具体链表结构
//                  图见潮流计算文档
// 变量说明：       1、linkcnt     为经验公式，意在开辟较大的链表空间
//                  2、colstatept  指针指向当前节点所连当前支路的尾节点   
// 注意事项：       链表中各链节众多指针的指向比较复杂，容易混淆，需注意
// 注释信息：       mdn  8.16
//************************************************************************
void EmsMatrix::opinitialize(struct opworkst *workpt)
{
   int i, j,linkcnt;
   unsigned int size; 
   struct opstatest *statept, *colstatept;
   struct oporderst *orderpt;    
   struct oplinkst *linkpt, *next_linkpt;
   struct mastatest *mastatept;
   struct maamtrxst *maamtrxpt;
 
   //add by qwbin 2006-2-11
   nmastate = mastate_vt.size();
   nmaamtrx = maamtrx_vt.size();    
   nmagmtrx = magmtrx_vt.size();    
   //nmagmtrx = 12*nmaamtrx;

   workpt->off_diagonalpt = nmagmtrx;
    
   size   = nmastate * sizeof( struct opstatest );
   workpt->statept = (struct opstatest *)malloc(size);
   memset ( (char *) workpt->statept, 0, size );

   statept  = workpt->statept;

   for (i = 1; i <= nmastate; i++, statept++)
   {
      mastatept = mastate_vt[i-1];
	  statept->stateitm = i;   //statept 重新编号 
      statept->maamtrx_itm = mastatept->maamtrx_itm;
      statept->maamtrx_cnt = mastatept->maamtrx_cnt;
      if ( mastatept->stat.active )
      {
         statept->stat.active = TRUE;
         //对角元素为0的标志
         if ( mastatept->diag == 0.0 )
   	 {
            statept->stat.zero_diag = TRUE;
            statept->stat.reorder   = TRUE;
            workpt->num_zero_diagonal++;
   	 }
      }
   }
  
   size   = nmastate * sizeof( struct oporderst );
   workpt->orderpt = (struct oporderst *)malloc( size );
   memset ( (char *) workpt->orderpt, 0, size );
   
   for (i = 1, orderpt = workpt->orderpt; i <= nmastate; i++, orderpt++ )
      orderpt->orderitm = i;
   
   workpt->last_orderpt = orderpt - 1;
   workpt->skip_orderpt = orderpt - AS_ORDER_SKIP;
   
   //我理解这里是个经验值,没有精确的计算意义,对于列向量多的矩阵,参数必须增加
   linkcnt = AS_ORDER_STATE_MULT*nmastate + AS_ORDER_AMTRX_MULT*nmaamtrx;
   size = linkcnt * sizeof ( struct oplinkst );

   workpt->linkpt = (struct oplinkst *)malloc( size );
   workpt->garbagept = workpt->linkpt;
   linkpt = workpt->linkpt;
   
   //链表形成
   for ( i = 0; i < linkcnt; i++ )
   {
      next_linkpt = linkpt + 1;
      linkpt->linkpt = next_linkpt;
      linkpt = next_linkpt;
   }
   linkpt--;

   //linkpt指向最后一个地址
   linkpt->linkpt = NULL;
   
   statept = workpt->statept;
  
   //以下代码的功能是:通过全部矩阵元素的遍历,找出与各节点的连接全部节点 
   // //对行循环statept++,改变了链表关系,指向前一LINKPT
   for ( i = 0; i < nmastate; i++, statept++ )  //按节点循环  mdn  8.2
   {
      if ( statept->stat.active )
      {
      //行首元素.  
	  //maamtrxpt = workpt->maamtrxpt + statept->maamtrx_itm - 1;
	 //对该行的列循环,形成linkpt链表(有元素存在,表示行,列标号的两节点关联)
     for ( j = 0; j < statept->maamtrx_cnt; j++)//按节点所连支路循环  mdn  8.2
   	 {
            maamtrxpt =maamtrx_vt[j+statept->maamtrx_itm-1];
		    //矩阵元素对应列节点指针
		    colstatept = workpt->statept + maamtrxpt->colitm - 1;//指向当前节点所连支路的末端节点！mdn  8.6
            if ( colstatept->stat.active )
            {
               //临时数据,前次linkpt记录
			   linkpt = workpt->garbagept; //临时变量，完成linkst结构链表的遍历  mdn  8.6
			   //记录下次linkpt值,下一次用
               workpt->garbagept = linkpt->linkpt;//已经构造完毕,由此实现linkpt遍列
               
			   //linkpt->statept  = 找到的矩阵元素对应列节点.
			   linkpt->statept = colstatept; //指向当前支路的尾节点  mdn  8.6 
               
			   //行的下一LINK值与STATE的现在LINK一致 //内循环中statept指针是不变的
			   linkpt->linkpt  = statept->linkpt; //形成前向链表  mdn  8.6
               
			   //返回给statept,记录下LINK指针
			   statept->linkpt = linkpt;//好象是为了变换指针,由此实现linkpt遍列,做临时变量用
               
			   //计数
			   statept->connection_cnt++;//当前节点的度  mdn  8.6
            }
         }
      }
   }
   return;
}

//************************************************************************
// test_symmetry - test for symmetric matrix
// 统计mastate->cnt 比较statept的linkpt数是否一致. 
//************************************************************************
// 函数名称：       test_symmetry 
// 返回值:          int
// 函数功能：       用来测试形成的链表结构是否成称，即链表是否与原有节点导纳阵
//                 （或网络）结构相同。
// 判断原则：       比较statept->linkpt（指向行链表首链节）所连链表的链节数目是
//                  否与statept->connection_cnt（节点导纳阵形成过程中计数）相等。
//                  如二者相等，则说明链表结构对称；不等，则说明链表结构不对称，
//                  需调用函数make_symmetric( )，以使链表结构对称。
// 变量说明：       
// 注意事项：       链表中各链节众多指针的指向比较复杂，容易混淆，需注意
// 注释信息：       mdn  8.16
//************************************************************************
int EmsMatrix::test_symmetry (struct opworkst *workpt)
{
   int i, nrow = 0, nhole = 0;
   int rc = SUCCESS;
   struct opstatest *statept;
   struct oplinkst *linkpt;
   
   statept = workpt->statept;
   for ( i = 0; i < nmastate; i++, statept++ )
   {
      linkpt = statept->linkpt;  //指向当前节点所连最后一条支路的链节  mdn  8.10
      while ( linkpt != NULL )  
      {
         linkpt->statept->cnt++;  //当前节点度加1  mdn  8.10
         linkpt  = linkpt->linkpt;  //移到下一链节，此处为内存地址小的链节即上一支路  mdn  8.10
      }
   }
   
   statept = workpt->statept;
   for ( i = 0; i < nmastate; i++, statept++ )
   {
      if ( statept->cnt != statept->connection_cnt )  //不对称，即链表形成过程有误  mdn  8.10
      {
         nrow++;
         rc = !SUCCESS;
         nhole += abs( statept->cnt - statept->connection_cnt );
      }
   }
   return ( rc );
}

//************************************************************************
// make_symmetric - establish symmetric matrix
// 形成对称矩阵 
//************************************************************************
void EmsMatrix::make_symmetric( struct opworkst *workpt)
{
   int i, nrow = 0, nhole = 0;
   struct opstatest *statept, *jstatept;
   struct oplinkst *linkpt, *jlinkpt; 
   
   statept = workpt->statept;
   
   //对节点循环
   for ( i = 0; i<nmastate; i++,statept++ )
   {
      linkpt = statept->linkpt;
      while ( linkpt != NULL )
      {
         jstatept = linkpt->statept;
         jlinkpt  = jstatept->linkpt;
         while ( jlinkpt != NULL )
         {
            if ( jlinkpt->statept == statept )
               break;
            jlinkpt = jlinkpt->linkpt;
         }
         linkpt = linkpt->linkpt;
   
         if ( jlinkpt == NULL )
         {
            if ( workpt->garbagept == NULL )
            {
               linkpt = NULL;
               workpt->stat.overflow = TRUE;
            }
            else
            {
               jstatept->stat.hole = TRUE;
               nhole++;
               jlinkpt = workpt->garbagept;
               workpt->garbagept = jlinkpt->linkpt;
               jlinkpt->statept = statept;
               jlinkpt->linkpt   = jstatept->linkpt;
               jstatept->linkpt  = jlinkpt;
            }
         }
      }
      if ( workpt->stat.overflow )
         break;
   }
   
   return;
}

//************************************************************************
// scheme1 - perform scheme I ordering
// 形成序，按自然顺序，形成orderpt
//************************************************************************
//************************************************************************
// 函数名称：       scheme1 
// 返回值:          无
// 函数功能：       为节点优化编号的主要步骤之一，按各节点度由小到大的顺序排列，
//                  通过numconpt指针（由base_numconpt处开始）找到度最小的节点链
//                  表linkst，再通过linkpt->statept找到度最小的节点。
// 变量说明：          
// 注意事项：       此处linkst链节含义（主要存储节点信息）不同于初始化中linkst链节
//                  含义（主要存储节点所连支路信息）。
// 注释信息：       mdn  8.16
//************************************************************************
void EmsMatrix::scheme1 ( struct opworkst *workpt )
{
   int i, max_numcon,norder, cnt;                     
   unsigned int size;         
   struct numconst *base_numconpt, *numconpt;     
   struct opstatest *statept;      
   struct oporderst *orderpt;      
   struct oplinkst *linkpt, *next_linkpt;  
   
   max_numcon = 0;
   //找出各节点的度的最大值
   for (i = 0, statept = workpt->statept; i < nmastate; i++, statept++ )//找出各节点的度的最大值 mdn 8.11
   {
      if ( statept->connection_cnt > max_numcon )
         max_numcon = statept->connection_cnt;
   }
   
   max_numcon++;
   size   = ( max_numcon + 1 ) * sizeof ( struct numconst );
   base_numconpt = (struct numconst *)malloc ( size );
   memset ( (char *) base_numconpt, 0, size );

   //对全部state循环,形成 numconpt，numconpt存放linkpt
   for (i = 0, statept = workpt->statept; i < nmastate; i++, statept++ )//形成的链表结构详见图1  mdn 8.11
   {
      linkpt = workpt->garbagept;
      workpt->garbagept = linkpt->linkpt;
      if ( statept->stat.active )
         cnt = statept->connection_cnt + 1;
      else
         cnt = 0;
      numconpt = base_numconpt + cnt;
       
      linkpt->statept  = statept;
      linkpt->linkpt   = numconpt->linkpt;
      numconpt->linkpt = linkpt;
   }
   
   orderpt = workpt->orderpt;
   norder = 1;
   //按自然顺序，形成orderpt，
   for (i = 0, numconpt = base_numconpt; i <= max_numcon; i++, numconpt++ )////形成的链表结构详见图1  mdn 8.11
   {
      linkpt = numconpt->linkpt;
      while ( linkpt != NULL )
      {
         statept = linkpt->statept;
         orderpt->statept = statept;
         statept->orderpt = orderpt;
         next_linkpt = linkpt->linkpt;
         linkpt->linkpt = workpt->garbagept;
         workpt->garbagept = linkpt;
         linkpt = next_linkpt;
         orderpt++; // 为节点优化后的链表，orderpt->statept->stateitm为优化后的顺序所对应的原始节点号  mdn  8.11
         norder++;  //相当于I5,即优化后的排序   mdn  8.10
      }
   }
   
   free ( base_numconpt );

   return;
}

//************************************************************************
// reorder_scheme1 - reorder Scheme I for zero diagonals
// 0对角元放到最前 
//************************************************************************
void EmsMatrix::reorder_scheme1 ( struct opworkst *workpt )
{
   int i;
 
   struct oporderst *orderpt, *new_orderpt, *last_orderpt;

   struct oplinkst *linkpt;
   
   //最后一行记录的指针
   last_orderpt = workpt->orderpt + nmastate - 1;
   
   orderpt = workpt->orderpt;
   for ( i = 0; i < nmastate; i++, orderpt++ )
   {
      if ( orderpt->statept->stat.active )
      {  
		 //0对角元 
         if ( orderpt->statept->stat.reorder )
         {
            new_orderpt = last_orderpt;
            linkpt = orderpt->statept->linkpt;
            while ( linkpt != NULL )
            {
       	      //链接不为0对角元 
			  if ( ! linkpt->statept->stat.reorder )
               {