src/img/imgMlp.c

Go to the documentation of this file.

#include "imgInt.h"
#include "fsm.h"

static char rcsid[] UNUSED = "$Id: imgMlp.c,v 1.31 2005/04/26 19:08:33 jinh Exp $";


/*---------------------------------------------------------------------------*/
/* Constant declarations                                                     */
/*---------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------*/
/* Type declarations                                                         */
/*---------------------------------------------------------------------------*/


/*---------------------------------------------------------------------------*/
/* Structure declarations                                                    */
/*---------------------------------------------------------------------------*/

typedef struct RcInfo_s {
  int   pos;    /* level */
  int   lNum;
  int   lMin;   /* index */
  int   lMax;   /* index */
  int   gNum;
  int   gMin;   /* index */
  int   gMax;   /* index */
  int   prevG;  /* level */
  int   nextG;  /* level */
  union {
    struct {
      bdd_t     *func;
      array_t   *nsVarBddArray;
    } row;
    struct {
      bdd_t     *var;
      int       type;
    } col;
  } data;
} RcInfo_t;

typedef struct RcList_s {
  int                   index;
  int                   otherIndex;
  struct RcList_s       *next;
  struct RcList_s       *prev;
} RcList_t;

typedef struct RcListInfo_s {
  int           num;
  int           maxNum;
  int           curIndex;
  st_table      *table;
  RcList_t      *cur;
  RcList_t      *head;
  RcList_t      *tail;
} RcListInfo_t;

typedef struct ClusterList_s {
  int                   flag;   /* 0 : to be conjuncted
                                   1 : start element
                                   2 : last element
                                   3 : isolated
                                */
  int                   start;
  int                   end;
  char                  *supports;
  int                   minCol; /* index */
  int                   maxCol; /* index */
  int                   nSupports;
  float                 affinity;
  int                   nQuantifyVars;
  mdd_t                 *product;
  array_t               *nsVarBddArray;
  struct ClusterList_s  *prev;
  struct ClusterList_s  *next;
} ClusterList_t;

typedef struct ClusterSortedList_s {
  ClusterList_t                 *list;
  struct ClusterSortedList_s    *next;
} ClusterSortedList_t;

typedef struct SccList_s {
  int                   startFunc;
  int                   lastFunc;
  int                   startVar;
  int                   lastVar;
  int                   nFuncs;
  int                   nVars;
  struct SccList_s      *next;
} SccList_t;

typedef struct LatchList_s {
  int           number;
  int           bddId;
  mdd_t         *relation;
  st_table      *table;
} LatchList_t;


/*---------------------------------------------------------------------------*/
/* Variable declarations                                                     */
/*---------------------------------------------------------------------------*/

static  int     nMoves;


/*---------------------------------------------------------------------------*/
/* Macro declarations                                                        */
/*---------------------------------------------------------------------------*/


/*---------------------------------------------------------------------------*/
/* Static function prototypes                                                */
/*---------------------------------------------------------------------------*/

static void SetupMlp(mdd_manager *mddManager, char **xy, int nRows, int nCols, int *rowOrder, int *colOrder, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *varPos, array_t *nsVarBddArray, int *nActiveRows, int *nActiveCols, array_t *nonAppearingVarBddArray, Img_DirectionType direction, ImgTrmOption_t *option);
static SccList_t * MlpDecomposeScc(mdd_manager *mddManager, char **xy, int nRows, int nActiveRows, int nActiveCols, int *rowOrder, int *colOrder, RcInfo_t *rowInfo, RcInfo_t *colInfo, int clusteredFlag, ImgTrmOption_t *option);
static int SccSortListIncreasingWithVars(const void *p1, const void *p2);
static int SccSortListDecreasingWithVars(const void *p1, const void *p2);
static int SccSortListIncreasingWithArea(const void *p1, const void *p2);
static int SccSortListDecreasingWithArea(const void *p1, const void *p2);
static int SccSortListIncreasingWithRatio(const void *p1, const void *p2);
static int SccSortListDecreasingWithRatio(const void *p1, const void *p2);
static int SccSortRc(const void *p1, const void *p2);
static void FreeSccList(SccList_t *sccHeadList);
static int NumOfSccs(SccList_t *sccHeadList);
static void BlockLowerTriangle(char **xy, int nRows, int nCols, int nActiveRows, int nActiveCols, SccList_t *sccList, int *rowOrder, int *colOrder, RcInfo_t *rowInfo, RcInfo_t *colInfo, ImgTrmOption_t *option);
static void MoveBestRows(char **xy, SccList_t *sccList, int nRows, int nCols, int *rowOrder, int *colOrder, RcInfo_t *rowInfo, RcInfo_t *colInfo, int startFunc, int lastFunc, int startVar, int lastVar, ImgTrmOption_t *option);
static void MoveBestCols(char **xy, SccList_t *sccList, int nRows, int nCols, int nActiveRows, int nActiveCols, int *rowOrder, int *colOrder, RcInfo_t *rowInfo, RcInfo_t *colInfo, int startFunc, int lastFunc, int startVar, int lastVar, ImgTrmOption_t *option);
static void MlpPostProcess(char **xy, SccList_t *sccList, int nVars, int nRows, int nCols, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder, Img_DirectionType direction, ImgTrmOption_t *option);
static float ComputeLambdaMlp(char **xy, int nVars, int nRows, int nCols, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder, Img_DirectionType direction, int mode, int asIs, ImgTrmOption_t *option);
static void FindAndMoveSingletonRows(char **xy, SccList_t *sccList, int nRows, int nCols, int *startFunc, int *lastFunc, int *startVar, int *lastVar, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder, ImgTrmOption_t *option);
static int MoveSingletonRow(char **xy, SccList_t *sccList, int nRows, int nCols, int x, int *startFunc, int *lastFunc, int *startVar, int *lastVar, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder, RcListInfo_t *candList, ImgTrmOption_t *option);
static void FindAndMoveSingletonCols(char **xy, SccList_t *sccList, int nRows, int nCols, int *startFunc, int *lastFunc, int *startVar, int *lastVar, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder, ImgTrmOption_t *option);
static int MoveSingletonCol(char **xy, SccList_t *sccList, int nRows, int nCols, int y, int *startFunc, int *lastFunc, int *startVar, int *lastVar, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder, RcListInfo_t *candList, ImgTrmOption_t *option);
static void MoveRowToTop(char **xy, int x, int startFunc, int lastFunc, int startVar, int lastVar, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder, RcListInfo_t *candList, int method);
static void MoveColToLeft(char **xy, int y, int startFunc, int lastFunc, int startVar, int lastVar, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder, RcListInfo_t *candList, int method);
static void MoveRowToBottom(char **xy, int x, int startFunc, int lastFunc, int startVar, int lastVar, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder, RcListInfo_t *candList, int method);
static void MoveColToRight(char **xy, int y, int startFunc, int lastFunc, int startVar, int lastVar, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder, RcListInfo_t *candList, int method);
static void PrintMatrix(char **xy, int nRows, int nCols, int *rowOrder, int *colOrder, RcInfo_t *rowInfo, RcInfo_t *colInfo, int startFunc, int lastFunc, int startVar, int lastVar);
static void PrintMatrixWithCluster(char **xy, ClusterList_t *headCluster, int nCols, int *rowOrder, int *colOrder, Img_DirectionType direction);
static void PrintClusterMatrix(ClusterList_t *headCluster, int nCols, int *colOrder, Img_DirectionType direction);
static void CheckMatrix(char **xy, SccList_t *sccList, int nRows, int nCols, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder, int startFunc, int lastFunc, int startVar, int lastVar, int local);
static void CheckCluster(ClusterList_t *headCluster, int nCols, RcInfo_t *colInfo, int *colOrder);
static void WriteMatrix(FILE *fout, char **xy, int nRows, int nCols, int *rowOrder, int *colOrder, RcInfo_t *rowInfo, RcInfo_t *colInfo);
static void PrintRow(char **xy, int nRows, int nCols, int *rowOrder, int *colOrder, int from, int to);
static void PrintCol(char **xy, int nRows, int nCols, int *rowOrder, int *colOrder, int from, int to);
static RcListInfo_t * SortedListAlloc(void);
static void SortedListFree(RcListInfo_t *candList);
static void SortedListInsert(RcListInfo_t *candList, int index, int otherIndex, RcInfo_t *otherInfo, int method);
static void SortedListDelete(RcListInfo_t *candList, int index);
static void SortedListMove(RcListInfo_t *candList, RcInfo_t *info, int index, int method);
static void CheckSortedList(RcListInfo_t *candList, RcInfo_t *info, int method);
static void MlpCluster(mdd_manager *mddManager, char **xy, int nRows, int nCols, int nActiveRows, int nActiveCols, int *nClusterRows, int *nClusterCols, int *rowOrder, int *colOrder, RcInfo_t *rowInfo, RcInfo_t *colInfo, array_t *clusterArray, array_t *arraySmoothVarBddArray, Img_DirectionType direction, int *cRowOrder, array_t *nsVarBddArray, int *sccBorder, int *varPos, ImgTrmOption_t *option);
static int MlpCountSupport(ClusterList_t *list, int *colOrder, int nActiveCols);
static float MlpSupportAffinity(ClusterList_t *curList, ClusterList_t *nextList, RcInfo_t *colInfo, int *colOrder, int nActiveCols, int clusterMethod);
static int RecursiveCluster(mdd_manager *mddManager, ClusterList_t *headCluster, ClusterSortedList_t *clusterSortedList, char **xy, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder, int nActiveRows, int nClusterCols, Img_DirectionType direction, int *varPos, int moveFlag, ImgTrmOption_t *option);
static int RemoveLocalVarsInCluster(mdd_manager *mddManager, char **xy, ClusterList_t *list, int nActiveRows, int nClusterCols, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder, int moveFlag, ImgTrmOption_t *option);
static int MlpNumQuantifyVars(ClusterList_t *list, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *colOrder, int nClusterCols);
static ClusterSortedList_t * ClusterSortedListInsert(ClusterSortedList_t *clusterSortedList, ClusterList_t *list, int useQuantifyVars);
static int CountClusterList(ClusterList_t *clusterList);
static int CountClusterSortedList(ClusterSortedList_t *clusterSortedList);
static array_t * CreateInitialCluster(mdd_manager *mddManager, array_t *relationArray, ImgFunctionData_t *functionData, array_t *nsVarBddArray, ImgTrmOption_t *option);
static void SortCol(char **xy, int nRows, int nCols, RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder);
static void UpdateDisapearingPsVars(mdd_manager *mddManager, char **xy, int nActiveRows, int nActiveCols, int *rowOrder, int *colOrder, RcInfo_t *rowInfo, RcInfo_t *colInfo, int row, ImgTrmOption_t *option);
static void UpdateDisapearingPsVarsInCluster(mdd_manager *mddManager, char **xy, int nActiveRows, int nActiveCols, int *rowOrder, int *colOrder, RcInfo_t *rowInfo, RcInfo_t *colInfo, ClusterList_t *list, int moveFlag, ImgTrmOption_t *option);
static void UpdateNonappearingNsVars(mdd_manager *mddManager, array_t *nsVarBddArray, int nRows, RcInfo_t *rowInfo, int *rowOrder, array_t *nonAppearingVarBddArray);
static void WriteOrder(FILE *fout, int nCols, int *colOrder, RcInfo_t *colInfo);

/*---------------------------------------------------------------------------*/
/* Definition of exported functions                                          */
/*---------------------------------------------------------------------------*/


void
Img_PrintMlpOptions(void)
{
  ImgTrmOption_t        *options;
  char                  dummy[80];

  options = ImgGetTrmOptions();

  switch (options->mlpCluster) {
  case 0 :
    sprintf(dummy, "linear clustering");
    break;
  case 1 :
    sprintf(dummy, "affinity based tree clustering (default)");
    break;
  default :
    sprintf(dummy, "invalid");
    break;
  }
  fprintf(vis_stdout, "MLP: mlp_cluster = %d (%s)\n",
          options->mlpCluster, dummy);

  switch (options->mlpReorder) {
  case 0 :
    sprintf(dummy, "no reordering after clustering (default)");
    break;
  case 1 :
    sprintf(dummy, "reorder using MLP after clustering (inside)");
    break;
  case 2 :
    sprintf(dummy, "reorder using MLP after clustering (outside)");
    break;
  case 3 :
    sprintf(dummy, "reorder using IWLS95 after clustering");
    break;
  default :
    sprintf(dummy, "invalid");
    break;
  }
  fprintf(vis_stdout, "MLP: mlp_reorder = %d (%s)\n",
          options->mlpReorder, dummy);

  switch (options->mlpPreReorder) {
  case 0 :
    sprintf(dummy, "no reordering after clustering (default)");
    break;
  case 1 :
    sprintf(dummy, "reorder using MLP after clustering (inside)");
    break;
  case 2 :
    sprintf(dummy, "reorder using MLP after clustering (outside)");
    break;
  case 3 :
    sprintf(dummy, "reorder using IWLS95 after clustering");
    break;
  default :
    sprintf(dummy, "invalid");
    break;
  }
  fprintf(vis_stdout, "MLP: mlp_pre_reorder = %d (%s)\n",
          options->mlpPreReorder, dummy);

  switch (options->mlpPostProcess) {
  case 0 :
    sprintf(dummy, "no postprocessing (default)");
    break;
  case 1 :
    sprintf(dummy, "do postprocessing after ordering");
    break;
  case 2 :
    sprintf(dummy, "do postprocessing after clustering or reordering");
    break;
  case 3 :
    sprintf(dummy, "do both 1 and 2");
    break;
  default :
    sprintf(dummy, "invalid");
    break;
  }
  fprintf(vis_stdout, "MLP: mlp_postprocess = %d (%s)\n",
          options->mlpPostProcess, dummy);
}


/*---------------------------------------------------------------------------*/
/* Definition of internal functions                                          */
/*---------------------------------------------------------------------------*/


bdd_t *
ImgMlpMultiwayAndSmooth(mdd_manager *mddManager,
                        ImgFunctionData_t *functionData,
                        array_t *relationArray,
                        array_t *domainVarBddArray,
                        array_t *quantifyVarBddArray,
                        array_t *rangeVarBddArray,
                        Img_DirectionType direction,
                        ImgTrmOption_t *option)
{
  int           i, clusterSize;
  array_t       *clusteredRelationArray;
  mdd_t         *result, *relation, *tmp;

  if (direction == Img_Forward_c) {
    if (array_n(domainVarBddArray) == 0 && array_n(quantifyVarBddArray) == 0) {
      if (array_n(relationArray) == 1) {
        relation = array_fetch(mdd_t *, relationArray, 0);
        result = bdd_dup(relation);
      } else {
        result = array_fetch(mdd_t *, relationArray, 0);
        for (i = 1; i < array_n(relationArray); i++) {
          relation = array_fetch(mdd_t *, relationArray, i);
          tmp = result;
          result = bdd_and(tmp, relation, 1, 1);
          mdd_free(tmp);
        }
      }
      return(result);
    }
  }

  clusterSize = option->clusterSize;
  option->clusterSize = 1000000000;
  ImgMlpClusterRelationArray(mddManager, functionData, relationArray,
                domainVarBddArray, quantifyVarBddArray, rangeVarBddArray,
                direction, &clusteredRelationArray, NIL(array_t *), option);
  option->clusterSize = clusterSize;

  assert(array_n(clusteredRelationArray) == 1);
  result = array_fetch(mdd_t *, clusteredRelationArray, 0);
  array_free(clusteredRelationArray);

  return(result);
}


void
ImgMlpClusterRelationArray(mdd_manager *mddManager,
                        ImgFunctionData_t *functionData,
                        array_t *relationArray,
                        array_t *domainVarBddArray,
                        array_t *quantifyVarBddArray,
                        array_t *rangeVarBddArray,
                        Img_DirectionType direction,
                        array_t **clusteredRelationArrayPtr,
                        array_t **arraySmoothVarBddArrayPtr,
                        ImgTrmOption_t *option)
{
  int           i, j, x, y, nRows, nCols, nActiveRows, nActiveCols;
  int           *rowOrder, *colOrder, *cRowOrder;
  RcInfo_t      *rowInfo, *colInfo;
  char          **xy;
  int           nClusterRows, nClusterCols;
  array_t       *clusterArray;
  bdd_t         *cluster, *relation, *tempCluster, *var, *nsVar;
  int           row, col, s, t;
  array_t       *arraySmoothVarBddArray, *smoothVarBddArray;
  int           qVarPos;
  array_t       *nonAppearingVarBddArray;
  int           *varPos;
  array_t       *psVarBddArray, *nsVarBddArray;
  int           index, nVars, nc;
  long          initialTime, finalTime;
  array_t       *clusteredRelationArray;
  SccList_t     *sccHeadList, *sccList;
  int           *sccBorder;
  float         lambda1, lambda2, lambda3;
  FILE          *fout;

  if (option->mlpVerbosity)
    initialTime = util_cpu_time();
  else
    initialTime = 0;

  psVarBddArray = domainVarBddArray;
  nsVarBddArray = rangeVarBddArray;

  if (option->mlpInitialCluster && functionData) {
    clusteredRelationArray = CreateInitialCluster(mddManager, relationArray,
                                                  functionData, nsVarBddArray,
                                                  option);
  } else
    clusteredRelationArray = relationArray;

  nRows = array_n(clusteredRelationArray);
  if (direction == Img_Forward_c)
    nCols = array_n(domainVarBddArray) + array_n(quantifyVarBddArray);
  else
    nCols = array_n(rangeVarBddArray) + array_n(quantifyVarBddArray);

  xy = ALLOC(char *, nRows);
  for (i = 0; i < nRows; i++) {
    xy[i] = ALLOC(char, nCols);
    memset(xy[i], 0, sizeof(char) * nCols);
  }

  rowOrder = ALLOC(int, nRows);
  for (i = 0; i < nRows; i++)
    rowOrder[i] = i;
  colOrder = ALLOC(int, nCols);
  for (i = 0; i < nCols; i++)
    colOrder[i] = i;

  rowInfo = ALLOC(RcInfo_t, nRows);
  memset(rowInfo, 0, sizeof(RcInfo_t) * nRows);
  colInfo = ALLOC(RcInfo_t, nCols);
  memset(colInfo, 0, sizeof(RcInfo_t) * nCols);
  nVars = bdd_num_vars(mddManager);
  varPos = ALLOC(int, nVars);
  memset(varPos, 0, sizeof(int) * nVars);

  for (i = 0; i < nRows; i++) {
    relation = array_fetch(bdd_t *, clusteredRelationArray, i);
    rowInfo[i].data.row.func = bdd_dup(relation);
  }

  nc = 0;
  for (i = 0; i < array_n(psVarBddArray); i++) {
    var = array_fetch(bdd_t *, psVarBddArray, i);
    colInfo[nc].data.col.var = var;
    colInfo[nc].data.col.type = 1;
    index = (int)bdd_top_var_id(var);
    varPos[index] = nc;
    nc++;
  }
  for (i = 0; i < array_n(quantifyVarBddArray); i++) {
    var = array_fetch(bdd_t *, quantifyVarBddArray, i);
    colInfo[nc].data.col.var = var;
    colInfo[nc].data.col.type = 2;
    index = (int)bdd_top_var_id(var);
    varPos[index] = nc;
    nc++;
  }

  if (arraySmoothVarBddArrayPtr)
    nonAppearingVarBddArray = array_alloc(bdd_t *, 0);
  else
    nonAppearingVarBddArray = NIL(array_t);

  SetupMlp(mddManager, xy, nRows, nCols, rowOrder, colOrder,
           rowInfo, colInfo, varPos, nsVarBddArray, &nActiveRows, &nActiveCols,
           nonAppearingVarBddArray, direction, option);

  if (nActiveRows == 0) {
    clusterArray = array_alloc(mdd_t *, 0);
    cluster = bdd_one(mddManager);
    for (i = nActiveRows; i < nRows; i++) {
      row = rowOrder[i];
      relation = rowInfo[row].data.row.func;
      tempCluster = bdd_and(cluster, relation, 1, 1);
      bdd_free(cluster);
      cluster = tempCluster;
    }
    array_insert_last(bdd_t *, clusterArray, cluster);
    *clusteredRelationArrayPtr = clusterArray;

    if (arraySmoothVarBddArrayPtr) {
      arraySmoothVarBddArray = array_alloc(array_t *, 0);
      *arraySmoothVarBddArrayPtr = arraySmoothVarBddArray;
      array_insert_last(array_t *, arraySmoothVarBddArray,
                        nonAppearingVarBddArray);
      if (direction == Img_Forward_c) {
        smoothVarBddArray = array_alloc(array_t *, 0);
        array_insert_last(array_t *, arraySmoothVarBddArray, smoothVarBddArray);
      } else {
        smoothVarBddArray = array_alloc(array_t *, 0);
        for (i = 0; i < nRows; i++) {
          row = rowOrder[i];
          for (j = 0; j < array_n(rowInfo[row].data.row.nsVarBddArray); j++) {
            nsVar = array_fetch(bdd_t *, rowInfo[row].data.row.nsVarBddArray,
                                j);
            array_insert_last(bdd_t *, smoothVarBddArray, bdd_dup(nsVar));
          }
        }
        array_insert_last(array_t *, arraySmoothVarBddArray,
                          smoothVarBddArray);
      }
    }

    FREE(varPos);
    for (i = 0; i < nRows; i++) {
      if (xy[i])
        FREE(xy[i]);
    }
    FREE(xy);
    FREE(rowOrder);
    FREE(colOrder);
    for (i = 0; i < nRows; i++) {
      bdd_free(rowInfo[i].data.row.func);
      if (rowInfo[i].data.row.nsVarBddArray)
        array_free(rowInfo[i].data.row.nsVarBddArray);
    }
    FREE(rowInfo);
    FREE(colInfo);
    return;
  }
  sccHeadList = MlpDecomposeScc(mddManager, xy, nRows, nActiveRows, nActiveCols,
                                rowOrder, colOrder, rowInfo, colInfo,
                                0, option);
  if (option->mlpVerbosity >= 3 && option->mlpPrintMatrix) {
    PrintMatrix(xy, nActiveRows, nActiveCols, rowOrder, colOrder,
                rowInfo, colInfo, 0, nActiveRows - 1, 0, nActiveCols - 1);
  }
  sccList = sccHeadList;
  while (sccList) {
    BlockLowerTriangle(xy, nRows, nCols, nActiveRows, nActiveCols, sccList,
                        rowOrder, colOrder, rowInfo, colInfo, option);
    sccList = sccList->next;
  }
  if (option->mlpVerbosity >= 2) {
    lambda1 = ComputeLambdaMlp(xy, nCols, nActiveRows, nActiveCols,
                                rowInfo, colInfo, rowOrder, colOrder,
                                direction, 0, 0, option);
    lambda2 = ComputeLambdaMlp(xy, nCols, nActiveRows, nActiveCols,
                                rowInfo, colInfo, rowOrder, colOrder,
                                direction, 1, 0, option);
    lambda3 = ComputeLambdaMlp(xy, nCols, nActiveRows, nActiveCols,
                                rowInfo, colInfo, rowOrder, colOrder,
                                direction, 2, 0, option);
    fprintf(vis_stdout, "Lambda after MLP = %f (%f, %f)\n",
            lambda1, lambda2, lambda3);
  }

  if (option->mlpVerbosity) {
    finalTime = util_cpu_time();
    fprintf(vis_stdout, "time for MLP = %10g\n",
           (double)(finalTime - initialTime) / 1000.0);
  }
  if (option->mlpVerbosity >= 3 && option->mlpPrintMatrix) {
    PrintMatrix(xy, nActiveRows, nActiveCols, rowOrder, colOrder,
                rowInfo, colInfo, 0, nActiveRows - 1, 0, nActiveCols - 1);
  }

  if (option->mlpPostProcess == 1 || option->mlpPostProcess == 3) {
    for (x = 0; x < nActiveRows; x++) {
      row = rowOrder[x];
      rowInfo[row].lNum = rowInfo[row].gNum;
      rowInfo[row].lMin = rowInfo[row].gMin;
      rowInfo[row].lMax = rowInfo[row].gMax;
      rowInfo[row].prevG = -1;
      rowInfo[row].nextG = nActiveCols;
    }
    for (y = 0; y < nActiveCols; y++) {
      col = colOrder[y];
      colInfo[col].lNum = colInfo[col].gNum;
      colInfo[col].lMin = colInfo[col].gMin;
      colInfo[col].lMax = colInfo[col].gMax;
      colInfo[col].prevG = -1;
      colInfo[col].nextG = nActiveRows;
    }

    sccList = sccHeadList;
    while (sccList) {
      MlpPostProcess(xy, sccList, nCols, nActiveRows, nActiveCols,
                   rowInfo, colInfo, rowOrder, colOrder, direction, option);
      sccList = sccList->next;
    }
    if (option->mlpVerbosity >= 3 && option->mlpPrintMatrix) {
      PrintMatrix(xy, nActiveRows, nActiveCols, rowOrder, colOrder,
                  rowInfo, colInfo, 0, nActiveRows - 1, 0, nActiveCols - 1);
    }
  }

  if (option->mlpWriteOrder) {
    fout = fopen(option->mlpWriteOrder, "w");
    if (fout) {
      WriteOrder(fout, nActiveCols, colOrder, colInfo);
      fclose(fout);
    } else {
      fprintf(vis_stderr, "** img error: can't open file [%s]\n",
              option->mlpWriteOrder);
    }
  }

  clusterArray = array_alloc(bdd_t *, 0);
  if (arraySmoothVarBddArrayPtr) {
    arraySmoothVarBddArray = array_alloc(array_t *, 0);
    array_insert_last(array_t *, arraySmoothVarBddArray,
                        nonAppearingVarBddArray);
  } else
    arraySmoothVarBddArray = NIL(array_t);

  if ((direction == Img_Forward_c && option->mlpReorder) ||
      (direction == Img_Backward_c && option->mlpPreReorder) ||
      option->mlpPostProcess >= 2) {
    cRowOrder = ALLOC(int, nActiveRows);
    for (i = 0; i < nActiveCols; i++) {
      col = colOrder[i];
      colInfo[col].lNum = 0;
      colInfo[col].lMin = nActiveRows;
      colInfo[col].lMax = -1;
    }
  } else
    cRowOrder = NIL(int);

  if (option->mlpClusterScc) {
    sccBorder = ALLOC(int, nRows);
    memset(sccBorder, 0, sizeof(int) * nRows);
    sccList = sccHeadList;
    while (sccList) {
      sccBorder[sccList->startFunc] = 1;
      sccBorder[sccList->lastFunc] = 2;
      sccList = sccList->next;
    }
  } else
    sccBorder = NIL(int);

  MlpCluster(mddManager, xy, nRows, nCols, nActiveRows, nActiveCols,
                &nClusterRows, &nClusterCols,
                rowOrder, colOrder, rowInfo, colInfo,
                clusterArray, arraySmoothVarBddArray,
                direction, cRowOrder, nsVarBddArray, sccBorder, varPos, option);

  if (sccBorder)
    FREE(sccBorder);

  if ((direction == Img_Forward_c && option->mlpReorder) ||
      (direction == Img_Backward_c && option->mlpPreReorder) ||
      option->mlpPostProcess >= 2) {
    if (option->mlpDecomposeScc && NumOfSccs(sccHeadList) > 1) {
      FreeSccList(sccHeadList);
      sccHeadList = MlpDecomposeScc(mddManager, xy, nRows,
                                nClusterRows, nClusterCols,
                                cRowOrder, colOrder, rowInfo, colInfo,
                                1, option);
    } else {
      sccHeadList->startFunc = 0;
      sccHeadList->lastFunc = nClusterRows - 1;
      sccHeadList->startVar = 0;
      sccHeadList->lastVar = nClusterCols - 1;
      sccHeadList->nFuncs = nClusterRows;
      sccHeadList->nVars = nClusterCols;
    }
    if (option->mlpVerbosity >= 3 && option->mlpPrintMatrix) {
      PrintMatrix(xy, nClusterRows, nClusterCols, cRowOrder, colOrder,
                  rowInfo, colInfo, 0, nClusterRows - 1, 0, nClusterCols - 1);
    }

    if ((direction == Img_Forward_c && option->mlpReorder) ||
        (direction == Img_Backward_c && option->mlpPreReorder)) {
      sccList = sccHeadList;
      while (sccList) {
        BlockLowerTriangle(xy, nRows, nCols, nClusterRows, nClusterCols,
                           sccList, cRowOrder, colOrder, rowInfo, colInfo,
                           option);
        sccList = sccList->next;
      }
      if (option->mlpVerbosity >= 2) {
        lambda1 = ComputeLambdaMlp(xy, nCols, nClusterRows, nClusterCols,
                                   rowInfo, colInfo, rowOrder, colOrder,
                                   direction, 0, 0, option);
        lambda2 = ComputeLambdaMlp(xy, nCols, nClusterRows, nClusterCols,
                                   rowInfo, colInfo, rowOrder, colOrder,
                                   direction, 1, 0, option);
        lambda3 = ComputeLambdaMlp(xy, nCols, nClusterRows, nClusterCols,
                                   rowInfo, colInfo, rowOrder, colOrder,
                                   direction, 2, 0, option);
        fprintf(vis_stdout, "Lambda after MLP-cluster = %f (%f, %f)\n",
                lambda1, lambda2, lambda3);
      }

      if (option->mlpVerbosity >= 3 && option->mlpPrintMatrix) {
        PrintMatrix(xy, nClusterRows, nClusterCols, cRowOrder, colOrder,
                    rowInfo, colInfo, 0, nClusterRows - 1, 0, nClusterCols - 1);
      }
    }

    if (option->mlpPostProcess >= 2) {
      for (x = 0; x < nClusterRows; x++) {
        row = rowOrder[x];
        rowInfo[row].lNum = rowInfo[row].gNum;
        rowInfo[row].lMin = rowInfo[row].gMin;
        rowInfo[row].lMax = rowInfo[row].gMax;
        rowInfo[row].prevG = -1;
        rowInfo[row].nextG = nClusterCols;
      }
      for (y = 0; y < nClusterCols; y++) {
        col = colOrder[y];
        colInfo[col].lNum = colInfo[col].gNum;
        colInfo[col].lMin = colInfo[col].gMin;
        colInfo[col].lMax = colInfo[col].gMax;
        colInfo[col].prevG = -1;
        colInfo[col].nextG = nClusterRows;
      }

      sccList = sccHeadList;
      while (sccList) {
        MlpPostProcess(xy, sccList, nCols, nClusterRows, nClusterCols,
                     rowInfo, colInfo, cRowOrder, colOrder, direction, option);
        sccList = sccList->next;
      }
      if (option->mlpVerbosity >= 3 && option->mlpPrintMatrix) {
        PrintMatrix(xy, nClusterRows, nClusterCols, cRowOrder, colOrder,
                    rowInfo, colInfo, 0, nClusterRows - 1, 0, nClusterCols - 1);
      }
    }

    *clusteredRelationArrayPtr = clusterArray;

    if (direction == Img_Forward_c) {
      for (i = nClusterRows - 1; i >= 0; i--) {
        row = cRowOrder[i];
        array_insert_last(bdd_t *, clusterArray, rowInfo[row].data.row.func);
        rowInfo[row].data.row.func = NIL(bdd_t);
      }

      if (nRows > nActiveRows) {
        cluster = bdd_one(mddManager);
        for (i = nActiveRows; i < nRows; i++) {
          row = rowOrder[i];
          relation = rowInfo[row].data.row.func;
          tempCluster = bdd_and(cluster, relation, 1, 1);
          bdd_free(cluster);
          cluster = tempCluster;
        }
        array_insert_last(bdd_t *, clusterArray, cluster);
      }
    } else {
      if (nRows > nActiveRows) {
        UpdateNonappearingNsVars(mddManager, nsVarBddArray, nClusterRows,
                                 rowInfo, cRowOrder, nonAppearingVarBddArray);

        cluster = bdd_one(mddManager);
        for (i = nActiveRows; i < nRows; i++) {
          row = rowOrder[i];
          relation = rowInfo[row].data.row.func;
          tempCluster = bdd_and(cluster, relation, 1, 1);
          bdd_free(cluster);
          cluster = tempCluster;
        }
        array_insert_last(bdd_t *, clusterArray, cluster);
      }

      for (i = 0; i < nClusterRows; i++) {
        row = cRowOrder[i];
        array_insert_last(bdd_t *, clusterArray, rowInfo[row].data.row.func);
        rowInfo[row].data.row.func = NIL(bdd_t);
      }
    }

    if (arraySmoothVarBddArrayPtr) {
      *arraySmoothVarBddArrayPtr = arraySmoothVarBddArray;

      if (direction == Img_Forward_c) {
        if (nCols > nClusterCols) {
          for (i = nClusterCols; i < nCols; i++) {
            col = colOrder[i];
            array_insert_last(bdd_t *, nonAppearingVarBddArray,
                                bdd_dup(colInfo[col].data.col.var));
          }
        }
        qVarPos = nClusterCols - 1;
        if (qVarPos >= 0) {
          col = colOrder[qVarPos];
          while (colInfo[col].gNum == 0) {
            array_insert_last(bdd_t *, nonAppearingVarBddArray,
                                bdd_dup(colInfo[col].data.col.var));
            if (qVarPos == 0)
              break;
            qVarPos--;
            col = colOrder[qVarPos];
          }
        }
        for (i = nClusterRows - 1; i >= 0; i--) {
          row = cRowOrder[i];
          smoothVarBddArray = array_alloc(array_t *, 0);
          col = colOrder[qVarPos];
          while (rowInfo[colInfo[col].gMin].pos == i) {
            array_insert_last(bdd_t *, smoothVarBddArray,
                                bdd_dup(colInfo[col].data.col.var));
            if (qVarPos == 0)
              break;
            qVarPos--;
            col = colOrder[qVarPos];
          }
          array_insert_last(array_t *, arraySmoothVarBddArray,
                            smoothVarBddArray);
        }

        if (nRows > nActiveRows) {
          smoothVarBddArray = array_alloc(array_t *, 0);
          array_insert_last(array_t *, arraySmoothVarBddArray,
                            smoothVarBddArray);
        }
      } else {
        if (nRows == nActiveRows) {
          UpdateNonappearingNsVars(mddManager, nsVarBddArray, nClusterRows,
                                   rowInfo, cRowOrder, nonAppearingVarBddArray);
        } else {
          smoothVarBddArray = array_alloc(array_t *, 0);
          for (i = nActiveRows; i < nRows; i++) {
            row = rowOrder[i];
            for (j = 0; j < array_n(rowInfo[row].data.row.nsVarBddArray); j++) {
              nsVar = array_fetch(bdd_t *, rowInfo[row].data.row.nsVarBddArray,
                                  j);
              array_insert_last(bdd_t *, smoothVarBddArray, bdd_dup(nsVar));
            }
          }
          array_insert_last(array_t *, arraySmoothVarBddArray,
                            smoothVarBddArray);
        }
        for (i = 0; i < nClusterRows; i++) {
          row = cRowOrder[i];
          smoothVarBddArray = array_alloc(array_t *, 0);
          for (j = 0; j < array_n(rowInfo[row].data.row.nsVarBddArray); j++) {
            nsVar = array_fetch(bdd_t *, rowInfo[row].data.row.nsVarBddArray,
                                j);
            array_insert_last(bdd_t *, smoothVarBddArray, bdd_dup(nsVar));
          }

          s = colInfo[rowInfo[row].gMin].pos;
          t = colInfo[rowInfo[row].gMax].pos;
          for (j = s; j <= t; j++) {
            col = colOrder[j];
            if (colInfo[col].data.col.type == 2) {
              if (rowInfo[colInfo[col].gMax].pos == i) {
                array_insert_last(bdd_t *, smoothVarBddArray,
                                  bdd_dup(colInfo[col].data.col.var));
              }
            }
          }
          array_insert_last(array_t *, arraySmoothVarBddArray,
                            smoothVarBddArray);
        }
      }
    }

    FREE(cRowOrder);

    if (option->mlpVerbosity) {
      finalTime = util_cpu_time();
      fprintf(vis_stdout, "time for MLP-clustering-reorder = %10g\n",
              (double)(finalTime - initialTime) / 1000.0);
    }
  } else {
    *clusteredRelationArrayPtr = clusterArray;
    if (arraySmoothVarBddArrayPtr)
      *arraySmoothVarBddArrayPtr = arraySmoothVarBddArray;

    if (option->mlpVerbosity) {
      finalTime = util_cpu_time();
      fprintf(vis_stdout, "time for MLP-clustering = %10g\n",
              (double)(finalTime - initialTime) / 1000.0);
    }
  }

  FreeSccList(sccHeadList);

  FREE(varPos);
  for (i = 0; i < nRows; i++) {
    if (xy[i])
      FREE(xy[i]);
  }
  FREE(xy);
  FREE(rowOrder);
  FREE(colOrder);
  for (i = 0; i < nRows; i++) {
    if (rowInfo[i].data.row.func)
      bdd_free(rowInfo[i].data.row.func);
    if (rowInfo[i].data.row.nsVarBddArray)
      array_free(rowInfo[i].data.row.nsVarBddArray);
  }
  FREE(rowInfo);
  FREE(colInfo);
  if (option->mlpInitialCluster)
    mdd_array_free(clusteredRelationArray);
}


void
ImgMlpOrderRelationArray(mdd_manager *mddManager,
                        array_t *relationArray,
                        array_t *domainVarBddArray,
                        array_t *quantifyVarBddArray,
                        array_t *rangeVarBddArray,
                        Img_DirectionType direction,
                        array_t **orderedRelationArrayPtr,
                        array_t **arraySmoothVarBddArrayPtr,
                        ImgTrmOption_t *option)
{
  int           i, j, x, y, nRows, nCols, nActiveRows, nActiveCols;
  int           *rowOrder, *colOrder;
  RcInfo_t      *rowInfo, *colInfo;
  char          **xy;
  array_t       *orderedArray;
  bdd_t         *relation, *var, *nsVar;
  int           row, col, s, t;
  array_t       *arraySmoothVarBddArray, *smoothVarBddArray;
  int           qVarPos;
  array_t       *nonAppearingVarBddArray;
  int           *varPos;
  array_t       *psVarBddArray, *nsVarBddArray;
  int           index, nVars, nc;
  long          initialTime, finalTime;
  SccList_t     *sccHeadList, *sccList;
  float         lambda1, lambda2, lambda3;

  arraySmoothVarBddArray = NIL(array_t);
  nRows = array_n(relationArray);
  if (direction == Img_Forward_c)
    nCols = array_n(domainVarBddArray) + array_n(quantifyVarBddArray);
  else
    nCols = array_n(rangeVarBddArray) + array_n(quantifyVarBddArray);

  if (nCols == 0) {
    if (direction == Img_Forward_c) {
      orderedArray = mdd_array_duplicate(relationArray);
      *orderedRelationArrayPtr = orderedArray;
      if (arraySmoothVarBddArrayPtr) {
        arraySmoothVarBddArray = array_alloc(array_t *, 0);
        *arraySmoothVarBddArrayPtr = arraySmoothVarBddArray;
        for (i = 0; i <= nRows; i++) {
          smoothVarBddArray = array_alloc(array_t *, 0);
          array_insert_last(array_t *, arraySmoothVarBddArray,
                            smoothVarBddArray);
        }
      }
    } else {
      orderedArray = array_alloc(mdd_t *, 0);
      array_insert_last(mdd_t *, orderedArray, mdd_one(mddManager));
      *orderedRelationArrayPtr = orderedArray;
      if (arraySmoothVarBddArrayPtr) {
        arraySmoothVarBddArray = array_alloc(array_t *, 0);
        *arraySmoothVarBddArrayPtr = arraySmoothVarBddArray;
        for (i = 0; i <= 1; i++) {
          smoothVarBddArray = array_alloc(array_t *, 0);
          array_insert_last(array_t *, arraySmoothVarBddArray,
                            smoothVarBddArray);
        }
      }
    }
    return;
  }

  if (option->mlpVerbosity)
    initialTime = util_cpu_time();
  else
    initialTime = 0;

  xy = ALLOC(char *, nRows);
  for (i = 0; i < nRows; i++) {
    xy[i] = ALLOC(char, nCols);
    memset(xy[i], 0, sizeof(char) * nCols);
  }

  rowOrder = ALLOC(int, nRows);
  for (i = 0; i < nRows; i++)
    rowOrder[i] = i;
  colOrder = ALLOC(int, nCols);
  for (i = 0; i < nCols; i++)
    colOrder[i] = i;

  rowInfo = ALLOC(RcInfo_t, nRows);
  memset(rowInfo, 0, sizeof(RcInfo_t) * nRows);
  colInfo = ALLOC(RcInfo_t, nCols);
  memset(colInfo, 0, sizeof(RcInfo_t) * nCols);
  nVars = bdd_num_vars(mddManager);
  varPos = ALLOC(int, nVars);
  memset(varPos, 0, sizeof(int) * nVars);

  for (i = 0; i < nRows; i++) {
    relation = array_fetch(bdd_t *, relationArray, i);
    rowInfo[i].data.row.func = bdd_dup(relation);
  }

  psVarBddArray = domainVarBddArray;
  nsVarBddArray = rangeVarBddArray;
  nc = 0;
  for (i = 0; i < array_n(psVarBddArray); i++) {
    var = array_fetch(bdd_t *, psVarBddArray, i);
    colInfo[nc].data.col.var = var;
    colInfo[nc].data.col.type = 1;
    index = (int)bdd_top_var_id(var);
    varPos[index] = nc;
    nc++;
  }
  for (i = 0; i < array_n(quantifyVarBddArray); i++) {
    var = array_fetch(bdd_t *, quantifyVarBddArray, i);
    colInfo[nc].data.col.var = var;
    colInfo[nc].data.col.type = 2;
    index = (int)bdd_top_var_id(var);
    varPos[index] = nc;
    nc++;
  }

  if (arraySmoothVarBddArrayPtr)
    nonAppearingVarBddArray = array_alloc(bdd_t *, 0);
  else
    nonAppearingVarBddArray = NIL(array_t);
  SetupMlp(mddManager, xy, nRows, nCols, rowOrder, colOrder,
           rowInfo, colInfo, varPos, nsVarBddArray, &nActiveRows, &nActiveCols,
           nonAppearingVarBddArray, direction, option);
  if (nActiveRows == 0) {
    orderedArray = array_alloc(bdd_t *, 0);
    if (arraySmoothVarBddArrayPtr) {
      arraySmoothVarBddArray = array_alloc(array_t *, 0);
      array_insert_last(array_t *, arraySmoothVarBddArray,
                        nonAppearingVarBddArray);
    }
    for (i = 0; i < nRows; i++) {
      row = rowOrder[i];
      relation = rowInfo[row].data.row.func;
      array_insert_last(bdd_t *, orderedArray, mdd_dup(relation));
      if (arraySmoothVarBddArrayPtr) {
        if (direction == Img_Forward_c) {
          smoothVarBddArray = array_alloc(array_t *, 0);
          array_insert_last(array_t *, arraySmoothVarBddArray,
                            smoothVarBddArray);
        } else {
          smoothVarBddArray = rowInfo[row].data.row.nsVarBddArray;
          rowInfo[row].data.row.nsVarBddArray = NIL(array_t);
          array_insert_last(array_t *, arraySmoothVarBddArray,
                            smoothVarBddArray);
        }
      }
    }
    *orderedRelationArrayPtr = orderedArray;
    if (arraySmoothVarBddArrayPtr)
      *arraySmoothVarBddArrayPtr = arraySmoothVarBddArray;

    FREE(varPos);
    for (i = 0; i < nRows; i++) {
      if (xy[i])
        FREE(xy[i]);
    }
    FREE(xy);
    FREE(rowOrder);
    FREE(colOrder);
    for (i = 0; i < nRows; i++) {
      bdd_free(rowInfo[i].data.row.func);
      if (rowInfo[i].data.row.nsVarBddArray)
        array_free(rowInfo[i].data.row.nsVarBddArray);
    }
    FREE(rowInfo);
    FREE(colInfo);
    return;
  }
  sccHeadList = MlpDecomposeScc(mddManager, xy, nRows, nActiveRows, nActiveCols,
                                rowOrder, colOrder, rowInfo, colInfo,
                                0, option);
  if (option->mlpVerbosity >= 3 && option->mlpPrintMatrix) {
    PrintMatrix(xy, nActiveRows, nActiveCols, rowOrder, colOrder,
                rowInfo, colInfo, 0, nActiveRows - 1, 0, nActiveCols - 1);
  }
  sccList = sccHeadList;
  while (sccList) {
    BlockLowerTriangle(xy, nRows, nCols, nActiveRows, nActiveCols, sccList,
                        rowOrder, colOrder, rowInfo, colInfo, option);
    sccList = sccList->next;
  }
  if (option->mlpVerbosity >= 2) {
    lambda1 = ComputeLambdaMlp(xy, nCols, nActiveRows, nActiveCols,
                                rowInfo, colInfo, rowOrder, colOrder,
                                direction, 0, 0, option);
    lambda2 = ComputeLambdaMlp(xy, nCols, nActiveRows, nActiveCols,
                                rowInfo, colInfo, rowOrder, colOrder,
                                direction, 1, 0, option);
    lambda3 = ComputeLambdaMlp(xy, nCols, nActiveRows, nActiveCols,
                                rowInfo, colInfo, rowOrder, colOrder,
                                direction, 2, 0, option);
    fprintf(vis_stdout, "Lambda after MLP = %f (%f, %f)\n",
            lambda1, lambda2, lambda3);
  }

  if (option->mlpVerbosity) {
    finalTime = util_cpu_time();
    fprintf(vis_stdout, "time for MLP = %10g\n",
           (double)(finalTime - initialTime) / 1000.0);
  }
  if (option->mlpVerbosity >= 3 && option->mlpPrintMatrix) {
    PrintMatrix(xy, nActiveRows, nActiveCols, rowOrder, colOrder,
                rowInfo, colInfo, 0, nActiveRows - 1, 0, nActiveCols - 1);
  }

  if (option->mlpPostProcess) {
    for (x = 0; x < nRows; x++) {
      row = rowOrder[x];
      rowInfo[row].lNum = rowInfo[row].gNum;
      rowInfo[row].lMin = rowInfo[row].gMin;
      rowInfo[row].lMax = rowInfo[row].gMax;
      rowInfo[row].prevG = -1;
      rowInfo[row].nextG = nCols;
    }
    for (y = 0; y < nCols; y++) {
      col = colOrder[y];
      colInfo[col].lNum = colInfo[col].gNum;
      colInfo[col].lMin = colInfo[col].gMin;
      colInfo[col].lMax = colInfo[col].gMax;
      colInfo[col].prevG = -1;
      colInfo[col].nextG = nRows;
    }

    sccList = sccHeadList;
    while (sccList) {
      MlpPostProcess(xy, sccList, nCols, nActiveRows, nActiveCols,
                   rowInfo, colInfo, rowOrder, colOrder, direction, option);
      sccList = sccList->next;
    }
    if (option->mlpVerbosity >= 3 && option->mlpPrintMatrix) {
      PrintMatrix(xy, nActiveRows, nActiveCols, rowOrder, colOrder,
                  rowInfo, colInfo, 0, nActiveRows - 1, 0, nActiveCols - 1);
    }
  }

  orderedArray = array_alloc(bdd_t *, 0);
  if (arraySmoothVarBddArrayPtr) {
    arraySmoothVarBddArray = array_alloc(array_t *, 0);
    array_insert_last(array_t *, arraySmoothVarBddArray,
                        nonAppearingVarBddArray);
  }

  assert(nActiveCols > 0);
  if (direction == Img_Forward_c) {
    if (arraySmoothVarBddArrayPtr) {
      if (nCols > nActiveCols) {
        for (i = nActiveCols; i < nCols; i++) {
          col = colOrder[i];
          array_insert_last(bdd_t *, nonAppearingVarBddArray,
                            bdd_dup(colInfo[col].data.col.var));
        }
      }
      qVarPos = nActiveCols - 1;
      if (qVarPos >= 0) {
        col = colOrder[qVarPos];
        while (colInfo[col].gNum == 0) {
          array_insert_last(bdd_t *, nonAppearingVarBddArray,
                            bdd_dup(colInfo[col].data.col.var));
          if (qVarPos == 0)
            break;
          qVarPos--;
          col = colOrder[qVarPos];
        }
      }
    } else
      qVarPos = 0; /* to avoid warning */
    for (i = nActiveRows - 1; i >= 0; i--) {
      row = rowOrder[i];
      if (arraySmoothVarBddArrayPtr) {
        smoothVarBddArray = array_alloc(array_t *, 0);
        col = colOrder[qVarPos];
        while (rowInfo[colInfo[col].gMin].pos == i) {
          index = (int)bdd_top_var_id(colInfo[col].data.col.var);
          array_insert_last(bdd_t *, smoothVarBddArray,
                                bdd_dup(colInfo[col].data.col.var));
          if (qVarPos == 0)
            break;
          qVarPos--;
          col = colOrder[qVarPos];
        }
        array_insert_last(array_t *, arraySmoothVarBddArray, smoothVarBddArray);
      }
      relation = bdd_dup(rowInfo[row].data.row.func);
      array_insert_last(bdd_t *, orderedArray, relation);
    }

    if (nRows > nActiveRows) {
      for (i = nActiveRows; i < nRows; i++) {
        row = rowOrder[i];
        relation = rowInfo[row].data.row.func;
        array_insert_last(bdd_t *, orderedArray, mdd_dup(relation));
        if (arraySmoothVarBddArrayPtr) {
          smoothVarBddArray = array_alloc(array_t *, 0);
          array_insert_last(array_t *, arraySmoothVarBddArray,
                            smoothVarBddArray);
        }
      }
    }
  } else {
    if (arraySmoothVarBddArrayPtr) {
      UpdateNonappearingNsVars(mddManager, nsVarBddArray, nActiveRows,
                                rowInfo, rowOrder, nonAppearingVarBddArray);
    }
    if (nRows > nActiveRows) {
      for (i = nActiveRows; i < nRows; i++) {
        row = rowOrder[i];
        relation = rowInfo[row].data.row.func;
        array_insert_last(bdd_t *, orderedArray, mdd_dup(relation));
        if (arraySmoothVarBddArrayPtr) {
          smoothVarBddArray = rowInfo[row].data.row.nsVarBddArray;
          rowInfo[row].data.row.nsVarBddArray = NIL(array_t);
          array_insert_last(array_t *, arraySmoothVarBddArray,
                            smoothVarBddArray);
        }
      }
    }
    for (i = 0; i < nActiveRows; i++) {
      row = rowOrder[i];
      if (arraySmoothVarBddArrayPtr) {
        smoothVarBddArray = array_alloc(array_t *, 0);
        for (j = 0; j < array_n(rowInfo[row].data.row.nsVarBddArray); j++) {
          nsVar = array_fetch(bdd_t *, rowInfo[row].data.row.nsVarBddArray, j);
          array_insert_last(bdd_t *, smoothVarBddArray, bdd_dup(nsVar));
        }

        s = colInfo[rowInfo[row].gMin].pos;
        t = colInfo[rowInfo[row].gMax].pos;
        for (j = s; j <= t; j++) {
          col = colOrder[j];
          if (colInfo[col].data.col.type == 2) {
            if (rowInfo[colInfo[col].gMax].pos == i) {
              array_insert_last(bdd_t *, smoothVarBddArray,
                                bdd_dup(colInfo[col].data.col.var));
            }
          }
        }
        array_insert_last(array_t *, arraySmoothVarBddArray, smoothVarBddArray);
      }
      relation = bdd_dup(rowInfo[row].data.row.func);
      array_insert_last(bdd_t *, orderedArray, relation);
    }
  }

  *orderedRelationArrayPtr = orderedArray;
  if (arraySmoothVarBddArrayPtr)
    *arraySmoothVarBddArrayPtr = arraySmoothVarBddArray;

  FreeSccList(sccHeadList);

  FREE(varPos);
  for (i = 0; i < nRows; i++) {
    if (xy[i])
      FREE(xy[i]);
  }
  FREE(xy);
  FREE(rowOrder);
  FREE(colOrder);
  for (i = 0; i < nRows; i++) {
    bdd_free(rowInfo[i].data.row.func);
    if (rowInfo[i].data.row.nsVarBddArray)
      array_free(rowInfo[i].data.row.nsVarBddArray);
  }
  FREE(rowInfo);
  FREE(colInfo);
}


float
ImgMlpComputeLambda(mdd_manager *mddManager,
                    array_t *relationArray,
                    array_t *domainVarBddArray,
                    array_t *quantifyVarBddArray,
                    array_t *rangeVarBddArray,
                    Img_DirectionType direction,
                    int mode,
                    int asIs,
                    int *prevArea,
                    float *improvedLambda,
                    ImgTrmOption_t *option)
{
  int           i, nRows, nCols, nActiveRows, nActiveCols;
  int           *rowOrder, *colOrder;
  RcInfo_t      *rowInfo, *colInfo;
  char          **xy;
  bdd_t         *relation, *var;
  int           *varPos;
  array_t       *psVarBddArray, *nsVarBddArray;
  int           index, nVars, nc;
  SccList_t     *sccHeadList, *sccList;
  float         lambda;
  int           area;

  nRows = array_n(relationArray);
  if (direction == Img_Forward_c)
    nCols = array_n(domainVarBddArray) + array_n(quantifyVarBddArray);
  else
    nCols = array_n(rangeVarBddArray) + array_n(quantifyVarBddArray);

  if (nCols == 0) {
    if (improvedLambda) {
      *improvedLambda = 0.0;
      *prevArea = 0;
    }
    return(0.0);
  }

  xy = ALLOC(char *, nRows);
  for (i = 0; i < nRows; i++) {
    xy[i] = ALLOC(char, nCols);
    memset(xy[i], 0, sizeof(char) * nCols);
  }

  rowOrder = ALLOC(int, nRows);
  for (i = 0; i < nRows; i++)
    rowOrder[i] = i;
  colOrder = ALLOC(int, nCols);
  for (i = 0; i < nCols; i++)
    colOrder[i] = i;

  rowInfo = ALLOC(RcInfo_t, nRows);
  memset(rowInfo, 0, sizeof(RcInfo_t) * nRows);
  colInfo = ALLOC(RcInfo_t, nCols);
  memset(colInfo, 0, sizeof(RcInfo_t) * nCols);
  nVars = bdd_num_vars(mddManager);
  varPos = ALLOC(int, nVars);
  memset(varPos, 0, sizeof(int) * nVars);

  for (i = 0; i < nRows; i++) {
    relation = array_fetch(bdd_t *, relationArray, i);
    rowInfo[i].data.row.func = bdd_dup(relation);
  }

  psVarBddArray = domainVarBddArray;
  nsVarBddArray = rangeVarBddArray;
  nc = 0;
  for (i = 0; i < array_n(psVarBddArray); i++) {
    var = array_fetch(bdd_t *, psVarBddArray, i);
    colInfo[nc].data.col.var = var;
    colInfo[nc].data.col.type = 1;
    index = (int)bdd_top_var_id(var);
    varPos[index] = nc;
    nc++;
  }
  for (i = 0; i < array_n(quantifyVarBddArray); i++) {
    var = array_fetch(bdd_t *, quantifyVarBddArray, i);
    colInfo[nc].data.col.var = var;
    colInfo[nc].data.col.type = 2;
    index = (int)bdd_top_var_id(var);
    varPos[index] = nc;
    nc++;
  }

  SetupMlp(mddManager, xy, nRows, nCols, rowOrder, colOrder,
           rowInfo, colInfo, varPos, nsVarBddArray, &nActiveRows, &nActiveCols,
           NIL(array_t), direction, option);
  if (nActiveRows == 0) {
    FREE(varPos);
    for (i = 0; i < nRows; i++) {
      if (xy[i])
        FREE(xy[i]);
    }
    FREE(xy);
    FREE(rowOrder);
    FREE(colOrder);
    for (i = 0; i < nRows; i++) {
      bdd_free(rowInfo[i].data.row.func);
      if (rowInfo[i].data.row.nsVarBddArray)
        array_free(rowInfo[i].data.row.nsVarBddArray);
    }
    FREE(rowInfo);
    FREE(colInfo);

    if (improvedLambda) {
      *improvedLambda = 0.0;
      *prevArea = 0;
    }
    return(0.0);
  }
  sccHeadList = MlpDecomposeScc(mddManager, xy, nRows, nActiveRows, nActiveCols,
                                rowOrder, colOrder, rowInfo, colInfo,
                                0, option);
  sccList = sccHeadList;
  while (sccList) {
    BlockLowerTriangle(xy, nRows, nCols, nActiveRows, nActiveCols, sccList,
                        rowOrder, colOrder, rowInfo, colInfo, option);
    sccList = sccList->next;
  }

  lambda = ComputeLambdaMlp(xy, nCols, nActiveRows, nActiveCols,
                            rowInfo, colInfo, rowOrder, colOrder,
                            direction, mode, asIs, option);

  FreeSccList(sccHeadList);

  FREE(varPos);
  for (i = 0; i < nRows; i++) {
    if (xy[i])
      FREE(xy[i]);
  }
  FREE(xy);
  FREE(rowOrder);
  FREE(colOrder);
  for (i = 0; i < nRows; i++) {
    bdd_free(rowInfo[i].data.row.func);
    if (rowInfo[i].data.row.nsVarBddArray)
      array_free(rowInfo[i].data.row.nsVarBddArray);
  }
  FREE(rowInfo);
  FREE(colInfo);

  if (option->mlpLambdaMode == 0)
    area = nActiveRows * nActiveCols;
  else
    area = nActiveRows * nCols;
  if (improvedLambda) {
    if (*prevArea)
      *improvedLambda = lambda * (float)area / (*prevArea);
    else
      *improvedLambda = 0.0;
    *prevArea = area;
  }
  return(lambda);
}


void
ImgMlpGetQuantificationSchedule(mdd_manager *mddManager,
                        array_t *relationArray,
                        array_t *domainVarBddArray,
                        array_t *quantifyVarBddArray,
                        array_t *rangeVarBddArray,
                        array_t **clusteredRelationArrayPtr,
                        array_t **arraySmoothVarBddArrayPtr,
                        Img_DirectionType direction,
                        ImgTrmOption_t *option)
{
  int           i, j, nRows, nCols, nActiveRows, nActiveCols;
  int           *rowOrder, *colOrder;
  RcInfo_t      *rowInfo, *colInfo;
  char          **xy;
  bdd_t         *relation, *var, *nsVar;
  int           *varPos;
  array_t       *psVarBddArray, *nsVarBddArray;
  int           index, nVars, nc;
  array_t       *arraySmoothVarBddArray, *smoothVarBddArray;
  int           qVarPos, row, col, s, t;
  array_t       *nonAppearingVarBddArray;
  array_t       *clusteredRelationArray;

  nRows = array_n(relationArray);
  if (direction == Img_Forward_c)
    nCols = array_n(domainVarBddArray) + array_n(quantifyVarBddArray);
  else
    nCols = array_n(rangeVarBddArray) + array_n(quantifyVarBddArray);

  xy = ALLOC(char *, nRows);
  for (i = 0; i < nRows; i++) {
    xy[i] = ALLOC(char, nCols);
    memset(xy[i], 0, sizeof(char) * nCols);
  }

  rowOrder = ALLOC(int, nRows);
  for (i = 0; i < nRows; i++)
    rowOrder[i] = i;
  colOrder = ALLOC(int, nCols);
  for (i = 0; i < nCols; i++)
    colOrder[i] = i;

  rowInfo = ALLOC(RcInfo_t, nRows);
  memset(rowInfo, 0, sizeof(RcInfo_t) * nRows);
  colInfo = ALLOC(RcInfo_t, nCols);
  memset(colInfo, 0, sizeof(RcInfo_t) * nCols);
  nVars = bdd_num_vars(mddManager);
  varPos = ALLOC(int, nVars);
  memset(varPos, 0, sizeof(int) * nVars);

  psVarBddArray = domainVarBddArray;
  nsVarBddArray = rangeVarBddArray;
  if (direction == Img_Forward_c) {
    for (i = 0; i < nRows; i++) {
      relation = array_fetch(bdd_t *, relationArray, nRows - 1 - i);
      rowInfo[i].data.row.func = bdd_dup(relation);
    }
  } else {
    for (i = 0; i < nRows; i++) {
      relation = array_fetch(bdd_t *, relationArray, i);
      rowInfo[i].data.row.func = bdd_dup(relation);
    }
  }
  nc = 0;
  for (i = 0; i < array_n(psVarBddArray); i++) {
    var = array_fetch(bdd_t *, psVarBddArray, i);
    colInfo[nc].data.col.var = var;
    colInfo[nc].data.col.type = 1;
    index = (int)bdd_top_var_id(var);
    varPos[index] = nc;
    nc++;
  }
  for (i = 0; i < array_n(quantifyVarBddArray); i++) {
    var = array_fetch(bdd_t *, quantifyVarBddArray, i);
    colInfo[nc].data.col.var = var;
    colInfo[nc].data.col.type = 2;
    index = (int)bdd_top_var_id(var);
    varPos[index] = nc;
    nc++;
  }

  nonAppearingVarBddArray = array_alloc(bdd_t *, 0);
  SetupMlp(mddManager, xy, nRows, nCols, rowOrder, colOrder,
           rowInfo, colInfo, varPos, nsVarBddArray, &nActiveRows, &nActiveCols,
           nonAppearingVarBddArray, direction, option);
  SortCol(xy, nActiveRows, nActiveCols, rowInfo, colInfo, rowOrder, colOrder);

  clusteredRelationArray = array_alloc(mdd_t *, 0);
  arraySmoothVarBddArray = array_alloc(array_t *, 0);
  if (direction == Img_Forward_c) {
    if (nCols > nActiveCols) {
      for (i = nActiveCols; i < nCols; i++) {
        col = colOrder[i];
        array_insert_last(bdd_t *, nonAppearingVarBddArray,
                          bdd_dup(colInfo[col].data.col.var));
      }
    }
    qVarPos = nActiveCols - 1;
    if (qVarPos >= 0) {
      col = colOrder[qVarPos];
      while (colInfo[col].gNum == 0) {
        array_insert_last(bdd_t *, nonAppearingVarBddArray,
                          bdd_dup(colInfo[col].data.col.var));
        if (qVarPos == 0)
          break;
        qVarPos--;
        col = colOrder[qVarPos];
      }
    }
    array_insert_last(array_t *, arraySmoothVarBddArray,
                        nonAppearingVarBddArray);
    for (i = nActiveRows - 1; i >= 0; i--) {
      row = rowOrder[i];
      smoothVarBddArray = array_alloc(array_t *, 0);
      if (rowInfo[row].gNum > 0) {
        col = colOrder[qVarPos];
        while (rowInfo[colInfo[col].gMin].pos == i) {
          index = (int)bdd_top_var_id(colInfo[col].data.col.var);
          array_insert_last(bdd_t *, smoothVarBddArray,
                                bdd_dup(colInfo[col].data.col.var));
          if (qVarPos == 0)
            break;
          qVarPos--;
          col = colOrder[qVarPos];
        }
      }
      array_insert_last(array_t *, arraySmoothVarBddArray, smoothVarBddArray);
      relation = bdd_dup(rowInfo[row].data.row.func);
      array_insert_last(bdd_t *, clusteredRelationArray, relation);
    }
    if (nRows > nActiveRows) {
      for (i = nActiveRows; i < nRows; i++) {
        row = rowOrder[i];
        smoothVarBddArray = array_alloc(array_t *, 0);
        array_insert_last(array_t *, arraySmoothVarBddArray, smoothVarBddArray);
        relation = bdd_dup(rowInfo[row].data.row.func);
        array_insert_last(bdd_t *, clusteredRelationArray, relation);
      }
    }
  } else {
    array_insert_last(array_t *, arraySmoothVarBddArray,
                        nonAppearingVarBddArray);
    if (nRows > nActiveRows) {
      for (i = nActiveRows; i < nRows; i++) {
        row = rowOrder[i];
        smoothVarBddArray = rowInfo[row].data.row.nsVarBddArray;
        rowInfo[row].data.row.nsVarBddArray = NIL(array_t);
        array_insert_last(array_t *, arraySmoothVarBddArray, smoothVarBddArray);
        relation = bdd_dup(rowInfo[row].data.row.func);
        array_insert_last(bdd_t *, clusteredRelationArray, relation);
      }
    }
    for (i = 0; i < nActiveRows; i++) {
      row = rowOrder[i];
      smoothVarBddArray = array_alloc(array_t *, 0);
      for (j = 0; j < array_n(rowInfo[row].data.row.nsVarBddArray); j++) {
        nsVar = array_fetch(bdd_t *, rowInfo[row].data.row.nsVarBddArray, j);
        array_insert_last(bdd_t *, smoothVarBddArray, bdd_dup(nsVar));
      }
      if (rowInfo[row].gNum > 0) {
        s = colInfo[rowInfo[row].gMin].pos;
        t = colInfo[rowInfo[row].gMax].pos;
        for (j = s; j <= t; j++) {
          col = colOrder[j];
          if (colInfo[col].data.col.type == 2) {
            if (rowInfo[colInfo[col].gMax].pos == i) {
              array_insert_last(bdd_t *, smoothVarBddArray,
                                bdd_dup(colInfo[col].data.col.var));
            }
          }
        }
      }
      array_insert_last(array_t *, arraySmoothVarBddArray, smoothVarBddArray);
      relation = bdd_dup(rowInfo[row].data.row.func);
      array_insert_last(bdd_t *, clusteredRelationArray, relation);
    }
  }

  FREE(varPos);
  for (i = 0; i < nRows; i++) {
    if (xy[i])
      FREE(xy[i]);
  }
  FREE(xy);
  FREE(rowOrder);
  FREE(colOrder);
  for (i = 0; i < nRows; i++) {
    bdd_free(rowInfo[i].data.row.func);
    if (rowInfo[i].data.row.nsVarBddArray)
      array_free(rowInfo[i].data.row.nsVarBddArray);
  }
  FREE(rowInfo);
  FREE(colInfo);

  *clusteredRelationArrayPtr = clusteredRelationArray;
  *arraySmoothVarBddArrayPtr = arraySmoothVarBddArray;
}


void
ImgMlpPrintDependenceMatrix(mdd_manager *mddManager,
                        array_t *relationArray,
                        array_t *domainVarBddArray,
                        array_t *quantifyVarBddArray,
                        array_t *rangeVarBddArray,
                        Img_DirectionType direction,
                        int printFlag, FILE *fout,
                        ImgTrmOption_t *option)
{
  int           i, nRows, nCols, nActiveRows, nActiveCols;
  int           *rowOrder, *colOrder;
  RcInfo_t      *rowInfo, *colInfo;
  char          **xy;
  bdd_t         *relation, *var;
  int           *varPos;
  array_t       *psVarBddArray, *nsVarBddArray;
  int           index, nVars, nc;
  float         lambda1, lambda2, lambda3;

  nRows = array_n(relationArray);
  if (direction == Img_Forward_c)
    nCols = array_n(domainVarBddArray) + array_n(quantifyVarBddArray);
  else
    nCols = array_n(rangeVarBddArray) + array_n(quantifyVarBddArray);

  xy = ALLOC(char *, nRows);
  for (i = 0; i < nRows; i++) {
    xy[i] = ALLOC(char, nCols);
    memset(xy[i], 0, sizeof(char) * nCols);
  }

  rowOrder = ALLOC(int, nRows);
  for (i = 0; i < nRows; i++)
    rowOrder[i] = i;
  colOrder = ALLOC(int, nCols);
  for (i = 0; i < nCols; i++)
    colOrder[i] = i;

  rowInfo = ALLOC(RcInfo_t, nRows);
  memset(rowInfo, 0, sizeof(RcInfo_t) * nRows);
  colInfo = ALLOC(RcInfo_t, nCols);
  memset(colInfo, 0, sizeof(RcInfo_t) * nCols);
  nVars = bdd_num_vars(mddManager);
  varPos = ALLOC(int, nVars);
  memset(varPos, 0, sizeof(int) * nVars);

  psVarBddArray = domainVarBddArray;
  nsVarBddArray = rangeVarBddArray;
  if (direction == Img_Forward_c) {
    for (i = 0; i < nRows; i++) {
      relation = array_fetch(bdd_t *, relationArray, nRows - 1 - i);
      rowInfo[i].data.row.func = bdd_dup(relation);
    }
  } else {
    for (i = 0; i < nRows; i++) {
      relation = array_fetch(bdd_t *, relationArray, i);
      rowInfo[i].data.row.func = bdd_dup(relation);
    }
  }
  nc = 0;
  for (i = 0; i < array_n(psVarBddArray); i++) {
    var = array_fetch(bdd_t *, psVarBddArray, i);
    colInfo[nc].data.col.var = var;
    colInfo[nc].data.col.type = 1;
    index = (int)bdd_top_var_id(var);
    varPos[index] = nc;
    nc++;
  }
  for (i = 0; i < array_n(quantifyVarBddArray); i++) {
    var = array_fetch(bdd_t *, quantifyVarBddArray, i);
    colInfo[nc].data.col.var = var;
    colInfo[nc].data.col.type = 2;
    index = (int)bdd_top_var_id(var);
    varPos[index] = nc;
    nc++;
  }

  SetupMlp(mddManager, xy, nRows, nCols, rowOrder, colOrder,
           rowInfo, colInfo, varPos, nsVarBddArray, &nActiveRows, &nActiveCols,
           NIL(array_t), direction, option);
  SortCol(xy, nActiveRows, nActiveCols, rowInfo, colInfo, rowOrder, colOrder);

  if (printFlag) {
    PrintMatrix(xy, nActiveRows, nActiveCols, rowOrder, colOrder,
                rowInfo, colInfo, 0, nActiveRows - 1, 0, nActiveCols - 1);
  }
  if (fout) {
    WriteMatrix(fout, xy, nActiveRows, nActiveCols, rowOrder, colOrder,
                rowInfo, colInfo);
  }
  lambda1 = ComputeLambdaMlp(xy, nCols, nActiveRows, nActiveCols,
                                rowInfo, colInfo, rowOrder, colOrder,
                                direction, 0, 0, option);
  lambda2 = ComputeLambdaMlp(xy, nCols, nActiveRows, nActiveCols,
                                rowInfo, colInfo, rowOrder, colOrder,
                                direction, 1, 0, option);
  lambda3 = ComputeLambdaMlp(xy, nCols, nActiveRows, nActiveCols,
                                rowInfo, colInfo, rowOrder, colOrder,
                                direction, 2, 0, option);
  fprintf(vis_stdout, "Lambda = %f (%f, %f)\n", lambda1, lambda2, lambda3);

  FREE(varPos);
  for (i = 0; i < nRows; i++) {
    if (xy[i])
      FREE(xy[i]);
  }
  FREE(xy);
  FREE(rowOrder);
  FREE(colOrder);
  for (i = 0; i < nRows; i++) {
    bdd_free(rowInfo[i].data.row.func);
    if (rowInfo[i].data.row.nsVarBddArray)
      array_free(rowInfo[i].data.row.nsVarBddArray);
  }
  FREE(rowInfo);
  FREE(colInfo);
}


void
ImgMlpWriteClusterFile(FILE *fout,
                        mdd_manager *mddManager,
                        array_t *relationArray,
                        array_t *psVarBddArray,
                        array_t *nsVarBddArray)
{
  int           i, j, k, nVars;
  int           bddId, mddId, id, index;
  mdd_t         *relation;
  mdd_t         *psVar, *nsVar, **nsVars, **ns2ps;
  array_t       *supportBddIdArray, *bddIdArray;
  char          *name;
  int           count;

  nVars = bdd_num_vars(mddManager);
  nsVars = ALLOC(bdd_t *, nVars);
  memset(nsVars, 0, sizeof(bdd_t *) * nVars);
  ns2ps = ALLOC(bdd_t *, nVars);
  memset(ns2ps, 0, sizeof(bdd_t *) * nVars);

  for (i = 0; i < array_n(nsVarBddArray); i++) {
    psVar = array_fetch(bdd_t *, psVarBddArray, i);
    nsVar = array_fetch(bdd_t *, nsVarBddArray, i);
    index = (int)bdd_top_var_id(nsVar);
    nsVars[index] = nsVar;
    ns2ps[index] = psVar;
  }

  count = 0;
  for (i = 0; i < array_n(relationArray); i++) {
    relation = array_fetch(mdd_t *, relationArray, i);
    supportBddIdArray = mdd_get_bdd_support_ids(mddManager, relation);
    psVar = NULL;
    for (j = 0; j < array_n(supportBddIdArray); j++) {
      bddId = array_fetch(int, supportBddIdArray, j);
      psVar = ns2ps[bddId];
      if (psVar) {
        break;
      }
    }
    if (psVar) {
      if (count == 0)
        fprintf(fout, "CLUSTER[%d] {\n", count);
      else
        fprintf(fout, "\nCLUSTER[%d] {\n", count);
      count++;
    }
    for (j = 0; j < array_n(supportBddIdArray); j++) {
      bddId = array_fetch(int, supportBddIdArray, j);
      psVar = ns2ps[bddId];
      if (psVar) {
        name = mdd_read_var_name(psVar);
        nsVar = nsVars[bddId];
        mddId = mdd_read_mdd_id(nsVar);
        bddIdArray = mdd_id_to_bdd_id_array(mddManager, mddId);
        for (k = 0; k < array_n(bddIdArray); k++) {
          id = array_fetch(int, bddIdArray, k);
          if (id == bddId) {
            fprintf(fout, "%s %d\n", name, k);
            break;
          }
        }
        array_free(bddIdArray);
      }
    }
    array_free(supportBddIdArray);
    fprintf(fout, "}\n");
  }

  FREE(nsVars);
  FREE(ns2ps);
}




void
ImgMlpReadClusterFile(FILE *fin, mdd_manager *mddManager,
                        ImgFunctionData_t *functionData,
                        array_t *relationArray,
                        array_t *psVarBddArray,
                        array_t *nsVarBddArray,
                        array_t *quantifyVarBddArray,
                        Img_DirectionType direction,
                        array_t **clusteredRelationArrayPtr,
                        array_t **arraySmoothVarBddArrayPtr,
                        ImgTrmOption_t *option)
{
  array_t       *clusteredRelationArray, *newRelationArray;
  array_t       *arraySmoothVarBddArray;
  int           i, j, k, nVars;
  long          bddId;
  int           mddId, id, index;
  mdd_t         *relation, *cluster, *newCluster;
  mdd_t         *var, *psVar, *nsVar, **nsVars, **ns2ps;
  array_t       *supportBddIdArray, *bddIdArray;
  char          *name;
  LatchList_t   **latchList;
  st_table      *intermediateTable, *nameTable, *quantifyVarsTable, *idTable;
  int           nRelVars, relBddId, relMddId;
  mdd_t         *relVar;
  char          *resolved;
  char          line[512], nameArray[512];
  st_generator  *stGen1, *stGen2;

  nVars = bdd_num_vars(mddManager);
  nsVars = ALLOC(bdd_t *, nVars);
  memset(nsVars, 0, sizeof(bdd_t *) * nVars);
  ns2ps = ALLOC(bdd_t *, nVars);
  memset(ns2ps, 0, sizeof(bdd_t *) * nVars);
  latchList = ALLOC(LatchList_t *, nVars);
  memset(latchList, 0, sizeof(LatchList_t *) * nVars);
  resolved = ALLOC(char, nVars);
  memset(resolved, 0, sizeof(char) * nVars);

  quantifyVarsTable = st_init_table(st_numcmp, st_numhash);
  for (i = 0; i < array_n(functionData->quantifyBddVars); i++) {
    var = array_fetch(mdd_t *, functionData->quantifyBddVars, i);
    index = (int)bdd_top_var_id(var);
    st_insert(quantifyVarsTable, (char *)(long)index, NIL(char));
  }

  nameTable = st_init_table(strcmp, st_strhash);
  for (i = 0; i < array_n(nsVarBddArray); i++) {
    psVar = array_fetch(bdd_t *, psVarBddArray, i);
    nsVar = array_fetch(bdd_t *, nsVarBddArray, i);
    index = (int)bdd_top_var_id(nsVar);
    nsVars[index] = nsVar;
    ns2ps[index] = psVar;
    name = mdd_read_var_name(psVar);
    mddId = mdd_read_mdd_id(nsVar);
    st_insert(nameTable, name, (char *)(long)mddId);

    index = (int)bdd_top_var_id(psVar);
    st_insert(quantifyVarsTable, (char *)(long)index, NIL(char));
  }

  i = 0;
  st_foreach_item_int(nameTable, stGen1, &name, &mddId) {
    printf("[%d] name = %s mdd = %d\n", i, name, mddId);
    i++;
  }

  intermediateTable = st_init_table(st_ptrcmp, st_ptrhash);
  for (i = 0; i < array_n(relationArray); i++) {
    relation = array_fetch(mdd_t *, relationArray, i);
    supportBddIdArray = mdd_get_bdd_support_ids(mddManager, relation);
    nRelVars = 0;
    psVar = NIL(bdd_t);
    nsVar = NIL(bdd_t);
    relVar = NIL(bdd_t);
    relBddId = -1;
    relMddId = -1;
    idTable = st_init_table(st_numcmp, st_numhash);
    for (j = 0; j < array_n(supportBddIdArray); j++) {
      bddId = (long) array_fetch(int, supportBddIdArray, j);
      if (st_lookup(quantifyVarsTable, (char *)bddId, NIL(char *)))
        continue;
      psVar = ns2ps[bddId];
      if (psVar) {
        if (relVar) {
          bdd_free(relVar);
          relVar = NIL(bdd_t);
        }
        nsVar = nsVars[bddId];
        relMddId = mdd_read_mdd_id(nsVar);
        relBddId = (int) bddId;
        break;
      } else {
        if (nRelVars > 0) {
          if (relVar)
            bdd_free(relVar);
        }
        relVar = bdd_var_with_index(mddManager, (int) bddId);
        relMddId = mdd_read_mdd_id(relVar);
        relBddId = (int) bddId;
        name = mdd_read_var_name(relVar);
        st_insert(nameTable, name, (char *)(long)relMddId);
        st_insert(idTable, (char *)bddId, NIL(char));
        if (nRelVars >= 1) {
          bdd_free(relVar);
          relVar = NIL(bdd_t);
        }
        nRelVars++;
      }
    }
    array_free(supportBddIdArray);
    if (nsVar || relVar) {
      bddId = (int) relBddId;
      mddId = relMddId;

      if (relVar)
        resolved[bddId] = 1;

      bddIdArray = mdd_id_to_bdd_id_array(mddManager, mddId);
      for (k = 0; k < array_n(bddIdArray); k++) {
        id = array_fetch(int, bddIdArray, k);
        if (id == bddId) {
          if (latchList[mddId]) {
            if (latchList[mddId]->number == 1) {
              latchList[mddId]->table = st_init_table(st_numcmp, st_numhash);
              st_insert(latchList[mddId]->table,
                        (char *)(long)latchList[mddId]->bddId,
                        (char *)latchList[mddId]->relation);
            }
            st_insert(latchList[mddId]->table, (char *)bddId,
                      (char *)relation);
          } else {
            latchList[mddId] = ALLOC(LatchList_t, 1);
            latchList[mddId]->bddId = (int) bddId;
            latchList[mddId]->relation = relation;
            latchList[mddId]->number = 1;
            latchList[mddId]->table = NIL(st_table);
          }
          break;
        }
      }
      array_free(bddIdArray);
      st_free_table(idTable);
    } else {
      st_insert(intermediateTable, (char *)relation, (char *)idTable);
    }
  }

  i = 0;
  st_foreach_item_int(nameTable, stGen1, &name, &mddId) {
    printf("[%d] name = %s mdd = %d\n", i, name, mddId);
    i++;
  }

  while (intermediateTable->num_entries > 0) {
    st_foreach_item(intermediateTable, stGen1, &relation, &idTable) {
      st_foreach_item(idTable, stGen2, &bddId, NIL(char *)) {
        if (resolved[bddId])
          st_delete(idTable, &bddId, NIL(char *));
      }
      if (idTable->num_entries == 1) {
        st_foreach_item(idTable, stGen2, &bddId, NULL);
        relVar = bdd_var_with_index(mddManager, (int) bddId);
        mddId = mdd_read_mdd_id(relVar);
        mdd_free(relVar);

        st_delete(intermediateTable, &relation, NULL);

        bddIdArray = mdd_id_to_bdd_id_array(mddManager, mddId);
        for (k = 0; k < array_n(bddIdArray); k++) {
          id = array_fetch(int, bddIdArray, k);
          if (id == bddId) {
            if (latchList[mddId]) {
              if (latchList[mddId]->number == 1) {
                latchList[mddId]->table = st_init_table(st_numcmp, st_numhash);
                st_insert(latchList[mddId]->table,
                          (char *)(long)latchList[mddId]->bddId,
                          (char *)latchList[mddId]->relation);
              }
              st_insert(latchList[mddId]->table, (char *)bddId,
                        (char *)relation);
            } else {
              latchList[mddId] = ALLOC(LatchList_t, 1);
              latchList[mddId]->bddId = (int) bddId;
              latchList[mddId]->relation = relation;
              latchList[mddId]->number = 1;
              latchList[mddId]->table = NIL(st_table);
            }
            break;
          }
        }
        array_free(bddIdArray);
      }
    }
  }
  st_free_table(intermediateTable);

  i = 0;
  st_foreach_item_int(nameTable, stGen1, &name, &mddId) {
    printf("[%d] name = %s mdd = %d\n", i, name, mddId);
    i++;
  }

  cluster = NIL(mdd_t);
  clusteredRelationArray = array_alloc(mdd_t *, 0);
  while (fgets(line, 512, fin)) {
    if (line[0] == '\n' || line[0] == '#')
      continue;
    if (strncmp(line, "CLUSTER", 7) == 0)
      cluster = mdd_one(mddManager);
    else if (strncmp(line, "}", 1) == 0)
      array_insert_last(mdd_t *, clusteredRelationArray, cluster);
    else {
      sscanf(line, "%s %d", nameArray, &id);
      if (st_lookup_int(nameTable, nameArray, &mddId) == 0)
        assert(0);
      if (latchList[mddId]) {
        if (latchList[mddId]->table) {
          bddIdArray = mdd_id_to_bdd_id_array(mddManager, mddId);
          assert(id < array_n(bddIdArray));
          bddId = (long) array_fetch(int, bddIdArray, id);
          if (st_lookup(latchList[mddId]->table, (char *)bddId, &relation)) {
            latchList[mddId]->number--;
            if (latchList[mddId]->number == 0) {
              st_free_table(latchList[mddId]->table);
              FREE(latchList[mddId]);
              latchList[mddId] = NIL(LatchList_t);
            }
          } else
            assert(0);
        } else {
          relation = latchList[mddId]->relation;
          FREE(latchList[mddId]);
          latchList[mddId] = NIL(LatchList_t);
        }
        newCluster = bdd_and(cluster, relation, 1, 1);
        mdd_free(cluster);
        cluster = newCluster;
      } else
        assert(0);
    }
  }

  for (i = 0; i < nVars; i++) {
    if (latchList[i]) {
      if (latchList[i]->table) {
        st_foreach_item(latchList[i]->table, stGen1, &bddId, &relation) {
          cluster = bdd_dup(relation);
          array_insert_last(mdd_t *, clusteredRelationArray, cluster);
        }
        st_free_table(latchList[i]->table);
      } else {
        cluster = bdd_dup(latchList[i]->relation);
        array_insert_last(mdd_t *, clusteredRelationArray, cluster);
      }
      FREE(latchList[i]);
    }
  }

  FREE(nsVars);
  FREE(ns2ps);
  FREE(resolved);
  FREE(latchList);

  st_free_table(nameTable);
  st_free_table(quantifyVarsTable);

  if (arraySmoothVarBddArrayPtr) {
    ImgMlpGetQuantificationSchedule(mddManager,
                clusteredRelationArray, psVarBddArray, quantifyVarBddArray,
                nsVarBddArray, &newRelationArray, &arraySmoothVarBddArray,
                direction, option);
    mdd_array_free(clusteredRelationArray);
    clusteredRelationArray = newRelationArray;
  } else
    arraySmoothVarBddArray = NIL(array_t);

  *clusteredRelationArrayPtr = clusteredRelationArray;
  *arraySmoothVarBddArrayPtr = arraySmoothVarBddArray;
}


/*---------------------------------------------------------------------------*/
/* Definition of static functions                                            */
/*---------------------------------------------------------------------------*/

static void
SetupMlp(mdd_manager *mddManager,
         char **xy, int nRows, int nCols, int *rowOrder, int *colOrder,
         RcInfo_t *rowInfo, RcInfo_t *colInfo, int *varPos,
         array_t *nsVarBddArray, int *nActiveRows, int *nActiveCols,
         array_t *nonAppearingVarBddArray,
         Img_DirectionType direction, ImgTrmOption_t *option)
{
  int           i, j, x, y, s, t;
  int           row, col, otherCol, nr, nc, index;
  long          initialTime, finalTime;
  bdd_t         *relation, *nsVar;
  array_t       *supportArray, *localVarBddArray;
  int           support;
  array_t       **arrayLocalVarBddArray;
  float         lambda1, lambda2, lambda3;
  bdd_t         **nsVarArray;
  int           nVars;

  if (option->mlpVerbosity >= 2)
    initialTime = util_cpu_time();
  else
    initialTime = 0;

  nVars = bdd_num_vars(mddManager);
  nsVarArray = ALLOC(bdd_t *, nVars);
  memset(nsVarArray, 0, sizeof(bdd_t *) * nVars);
  for (i = 0; i < array_n(nsVarBddArray); i++) {
    nsVar = array_fetch(bdd_t *, nsVarBddArray, i);
    index = (int)bdd_top_var_id(nsVar);
    nsVarArray[index] = nsVar;
  }

  for (i = 0; i < nRows; i++) {
    relation = rowInfo[i].data.row.func;
    rowInfo[i].data.row.nsVarBddArray = array_alloc(bdd_t *, 0);
    supportArray = mdd_get_bdd_support_ids(mddManager, relation);
    for (j = 0; j < array_n(supportArray); j++) {
      support = array_fetch(int, supportArray, j);
      nsVar = nsVarArray[support];
      if (nsVar) {
        array_insert_last(bdd_t *, rowInfo[i].data.row.nsVarBddArray, nsVar);
        if (varPos[support])
          xy[i][varPos[support]] = 1;
        else {
          index = (int)bdd_top_var_id(colInfo[0].data.col.var);
          if (index == support)
            xy[i][varPos[support]] = 1;
        }
      } else
        xy[i][varPos[support]] = 1;
    }
    array_free(supportArray);
  }
  FREE(nsVarArray);

  for (x = 0; x < nRows; x++) {
    rowInfo[x].gMin = nCols;
    rowInfo[x].gMax = -1;
    rowInfo[x].pos = x;
    rowInfo[x].prevG = -1;
    rowInfo[x].nextG = nCols;
  }
  for (y = 0; y < nCols; y++) {
    colInfo[y].gMin = nRows;
    colInfo[y].gMax = -1;
    colInfo[y].pos = y;
    colInfo[y].prevG = -1;
    colInfo[y].nextG = nRows;
  }

  for (x = 0; x < nRows; x++) {
    for (y = 0; y < nCols; y++) {
      if (xy[x][y]) {
        if (x < colInfo[y].gMin)
          colInfo[y].gMin = x;
        if (x > colInfo[y].gMax)
          colInfo[y].gMax = x;
        if (y < rowInfo[x].gMin)
          rowInfo[x].gMin = y;
        if (y > rowInfo[x].gMax)
          rowInfo[x].gMax = y;
        rowInfo[x].gNum++;
        colInfo[y].gNum++;
      }
    }
  }

  nc = nCols;
  arrayLocalVarBddArray = ALLOC(array_t *, nRows);
  memset(arrayLocalVarBddArray, 0, sizeof(array_t *) * nRows);
  for (y = 0; y < nc; y++) {
    col = colOrder[y];
    if (colInfo[col].data.col.type == 2 && colInfo[col].gNum == 1) {
      row = colInfo[col].gMin;
      rowInfo[row].gNum--;
      if (rowInfo[row].gNum == 1) {
        if (rowInfo[row].gMin == col) {
          rowInfo[row].gMin = rowInfo[row].gMax;
          rowInfo[row].lMin = rowInfo[row].lMax;
        } else {
          rowInfo[row].gMax = rowInfo[row].gMin;
          rowInfo[row].lMax = rowInfo[row].lMin;
        }
      } else if (rowInfo[row].gNum > 1) {
        if (rowInfo[row].gMin == col) {
          s = colInfo[rowInfo[row].gMin].pos;
          t = colInfo[rowInfo[row].gMax].pos;
          for (j = s + 1; j <= t; j++) {
            otherCol = colOrder[j];
            if (xy[row][otherCol]) {
              rowInfo[row].gMin = otherCol;
              rowInfo[row].lMin = otherCol;
              break;
            }
          }
        } else if (rowInfo[row].gMax == col) {
          s = colInfo[rowInfo[row].gMin].pos;
          t = colInfo[rowInfo[row].gMax].pos;
          for (j = t - 1; j >= s; j--) {
            otherCol = colOrder[j];
            if (xy[row][otherCol]) {
              rowInfo[row].gMax = otherCol;
              rowInfo[row].lMax = otherCol;
              break;
            }
          }
        }
      }

      for (j = y; j < nc - 1; j++) {
        colOrder[j] = colOrder[j + 1];
        colInfo[colOrder[j]].pos = j;
      }
      colOrder[nc - 1] = col;
      colInfo[col].pos = nc - 1;

      nc--;
      y--;

      if (!arrayLocalVarBddArray[row])
        arrayLocalVarBddArray[row] = array_alloc(bdd_t *, 0);
      array_insert_last(bdd_t *, arrayLocalVarBddArray[row],
                        colInfo[col].data.col.var);
    }
  }

  if (direction == Img_Backward_c) {
    long        lindex;
    st_table    *unusedNsTable;
    st_generator *gen;

    unusedNsTable = st_init_table(st_numcmp, st_numhash);
    for (i = 0; i < array_n(nsVarBddArray); i++) {
      nsVar = array_fetch(bdd_t *, nsVarBddArray, i);
      lindex = (long) bdd_top_var_id(nsVar);
      st_insert(unusedNsTable, (char *)lindex, (char *)nsVar);
    }
    for (i = 0; i < nRows; i++) {
      for (j = 0; j < array_n(rowInfo[i].data.row.nsVarBddArray); j++) {
        nsVar = array_fetch(bdd_t *, rowInfo[i].data.row.nsVarBddArray, j);
        lindex = (long) bdd_top_var_id(nsVar);
        st_delete(unusedNsTable, &lindex, NULL);
      }
    }
    st_foreach_item(unusedNsTable, gen, &lindex, &nsVar) {
      array_insert_last(bdd_t *, nonAppearingVarBddArray, bdd_dup(nsVar));
    }
    st_free_table(unusedNsTable);
  }

  for (i = 0; i < nRows; i++) {
    localVarBddArray = arrayLocalVarBddArray[i];
    if (localVarBddArray) {
      relation = bdd_smooth(rowInfo[i].data.row.func, localVarBddArray);
      if (nonAppearingVarBddArray && direction == Img_Backward_c &&
          bdd_is_tautology(relation, 1)) {
          for (j = 0; j < array_n(rowInfo[i].data.row.nsVarBddArray); j++) {
            nsVar = array_fetch(bdd_t *, rowInfo[i].data.row.nsVarBddArray, j);
            array_insert_last(bdd_t *, nonAppearingVarBddArray, bdd_dup(nsVar));
          }
      }
      bdd_free(rowInfo[i].data.row.func);
      rowInfo[i].data.row.func = relation;
      UpdateDisapearingPsVars(mddManager, xy, nRows, nCols,
                        rowOrder, colOrder, rowInfo, colInfo,
                        i, option);
      if (localVarBddArray)
        array_free(localVarBddArray);
    }
  }
  FREE(arrayLocalVarBddArray);

  for (y = 0; y < nc; y++) {
    col = colOrder[y];
    if (colInfo[col].gNum == 0) {
      if (nonAppearingVarBddArray &&
          direction == Img_Forward_c && colInfo[col].data.col.type == 1) {
        array_insert_last(bdd_t *, nonAppearingVarBddArray,
                          bdd_dup(colInfo[col].data.col.var));
      }
      for (j = y; j < nc - 1; j++) {
        colOrder[j] = colOrder[j + 1];
        colInfo[colOrder[j]].pos = j;
      }
      colOrder[nc - 1] = col;
      colInfo[col].pos = nc - 1;
      nc--;
      y--;
    }
  }

  nr = nRows;
  for (x = 0; x < nr; x++) {
    row = rowOrder[x];
    if (rowInfo[row].gNum == 0) {
      for (i = x; i < nr - 1; i++) {
        rowOrder[i] = rowOrder[i + 1];
        rowInfo[rowOrder[i]].pos = i;
      }
      rowOrder[nr - 1] = row;
      rowInfo[row].pos = nr - 1;
      nr--;
      x--;
    }
  }

  for (x = 0; x < nr; x++) {
    row = rowOrder[x];
    rowInfo[row].lMin = rowInfo[row].gMin;
    rowInfo[row].lMax = rowInfo[row].gMax;
    rowInfo[row].lNum = rowInfo[row].gNum;
    if (nc != nCols && rowInfo[row].nextG == nCols)
      rowInfo[row].nextG = nc;
  }
  for (y = 0; y < nc; y++) {
    col = colOrder[y];
    colInfo[col].lMin = colInfo[col].gMin;
    colInfo[col].lMax = colInfo[col].gMax;
    colInfo[col].lNum = colInfo[col].gNum;
    if (nr != nRows && colInfo[col].nextG == nRows)
      colInfo[col].nextG = nr;
  }
  if (option->mlpVerbosity >= 2) {
    finalTime = util_cpu_time();
    fprintf(vis_stdout, "time for setup = %10g\n",
           (double)(finalTime - initialTime) / 1000.0);

    lambda1 = ComputeLambdaMlp(xy, nCols, nr, nc, rowInfo, colInfo,
                                rowOrder, colOrder, direction, 0, 0, option);
    lambda2 = ComputeLambdaMlp(xy, nCols, nr, nc, rowInfo, colInfo,
                                rowOrder, colOrder, direction, 1, 0, option);
    lambda3 = ComputeLambdaMlp(xy, nCols, nr, nc, rowInfo, colInfo,
                                rowOrder, colOrder, direction, 2, 0, option);
    fprintf(vis_stdout, "Initial Lambda = %f (%f, %f)\n",
            lambda1, lambda2, lambda3);
  }

  if (option->mlpDebug) {
    CheckMatrix(xy, NIL(SccList_t), nr, nc, rowInfo, colInfo,
                rowOrder, colOrder, 0, nr - 1, 0, nc - 1, 1);
  }

  *nActiveRows = nr;
  *nActiveCols = nc;
}


static SccList_t *
MlpDecomposeScc(mdd_manager *mddManager, char **xy, int nRows,
                int nActiveRows, int nActiveCols, int *rowOrder, int *colOrder,
                RcInfo_t *rowInfo, RcInfo_t *colInfo,
                int clusteredFlag, ImgTrmOption_t *option)
{
  int           i, j, x, y, size, nVars;
  int           startRow, startCol;
  int           nGroups, nChosen;
  int           row, col, otherRow;
  int           s1, t1, s2, t2;
  char          *rowFlag, *colFlag, *stackFlag;
  int           *stack;
  SccList_t     *sccHeadList, *sccTailList, *sccList;
  array_t       *rowScc, *colScc, *sccArray;
  long          initialTime, finalTime;
  int           *newRowOrder, *newColOrder;
  int           nCurRows, nCurCols;

  if (!option->mlpDecomposeScc) {
    sccList = ALLOC(SccList_t, 1);
    sccList->nFuncs = nActiveRows;
    sccList->nVars = nActiveCols;
    sccList->startFunc = 0;
    sccList->lastFunc = nActiveRows - 1;
    sccList->startVar = 0;
    sccList->lastVar = nActiveCols - 1;
    sccList->next = NIL(SccList_t);
    return(sccList);
  }

  if (option->mlpVerbosity >= 2)
    initialTime = util_cpu_time();
  else
    initialTime = 0;

  sccHeadList = NIL(SccList_t);
  sccTailList = NIL(SccList_t);

  rowFlag = ALLOC(char, nRows);
  memset(rowFlag, 0, sizeof(char) * nRows);
  nVars = bdd_num_vars(mddManager);
  colFlag = ALLOC(char, nVars);
  memset(colFlag, 0, sizeof(char) * nVars);
  stack = ALLOC(int, nActiveRows);
  stackFlag = ALLOC(char, nRows);
  memset(stackFlag, 0, sizeof(char) * nRows);

  startRow = 0;
  startCol = 0;
  nGroups = 0;
  nChosen = 0;
  while (nChosen < nActiveRows) {
    rowScc = array_alloc(int, 0);
    colScc = array_alloc(int, 0);

    size = 0;
    for (i = startRow; i < nActiveRows; i++) {
      row = rowOrder[i];
      if (rowFlag[row] == 0) {
        stack[0] = row;
        size = 1;
        stackFlag[row] = 1;
        break;
      }
    }

    while (size) {
      if (size + array_n(rowScc) == nRows) {
        FREE(stack);
        FREE(stackFlag);
        FREE(rowFlag);
        FREE(colFlag);
        array_free(rowScc);
        array_free(colScc);

        sccList = ALLOC(SccList_t, 1);
        sccList->nFuncs = nActiveRows;
        sccList->nVars = nActiveCols;
        sccList->startFunc = 0;
        sccList->lastFunc = nActiveRows - 1;
        sccList->startVar = 0;
        sccList->lastVar = nActiveCols - 1;
        sccList->next = NIL(SccList_t);
        return(sccList);
      }

      size--;
      row = stack[size];
      x = rowInfo[row].pos;
      rowFlag[row] = nGroups + 1;
      array_insert_last(int, rowScc, x);
      nChosen++;
      s1 = colInfo[rowInfo[row].gMin].pos;
      t1 = colInfo[rowInfo[row].gMax].pos;
      for (y = s1; y <= t1; y++) {
        col = colOrder[y];
        if (colInfo[col].gNum == nRows) {
          FREE(stack);
          FREE(stackFlag);
          FREE(rowFlag);
          FREE(colFlag);
          array_free(rowScc);
          array_free(colScc);

          sccList = ALLOC(SccList_t, 1);
          sccList->nFuncs = nActiveRows;
          sccList->nVars = nActiveCols;
          sccList->startFunc = 0;
          sccList->lastFunc = nActiveRows - 1;
          sccList->startVar = 0;
          sccList->lastVar = nActiveCols - 1;
          sccList->next = NIL(SccList_t);
          return(sccList);
        }
        if (!xy[row][col] || colFlag[col])
          continue;
        colFlag[col] = nGroups + 1;
        array_insert_last(int, colScc, y);
        s2 = rowInfo[colInfo[col].gMin].pos;
        t2 = rowInfo[colInfo[col].gMax].pos;
        for (i = s2; i <= t2; i++) {
          otherRow = rowOrder[i];
          if (rowFlag[otherRow] || stackFlag[otherRow])
            continue;
          if (xy[otherRow][col]) {
            stack[size] = otherRow;
            size++;
            stackFlag[otherRow] = 1;
          }
        }
      }
    }

    nGroups++;
    if (nGroups == 1 && nChosen == nActiveRows) {
      sccList = ALLOC(SccList_t, 1);
      sccList->nFuncs = nActiveRows;
      sccList->nVars = nActiveCols;
      sccList->startFunc = 0;
      sccList->lastFunc = nActiveRows - 1;
      sccList->startVar = 0;
      sccList->lastVar = nActiveCols - 1;
      sccList->next = NIL(SccList_t);

      sccHeadList = sccList;

      array_free(rowScc);
      array_free(colScc);
      break;
    }

    /* Move the rows and columns that belong to the current group to top-left */
    array_sort(rowScc, SccSortRc);
    array_sort(colScc, SccSortRc);
    for (i = 0; i < array_n(rowScc); i++) {
      x = array_fetch(int, rowScc, i);
      if (x == startRow + i)
        continue;
      row = rowOrder[x];
      for (j = x; j > startRow + i; j--) {
        rowOrder[j] = rowOrder[j - 1];
        rowInfo[rowOrder[j]].pos = j;
      }
      rowOrder[startRow + i] = row;
      rowInfo[row].pos = startRow + i;
    }
    for (i = 0; i < array_n(colScc); i++) {
      y = array_fetch(int, colScc, i);
      if (y == startCol + i)
        continue;
      col = colOrder[y];
      for (j = y; j > startCol + i; j--) {
        colOrder[j] = colOrder[j - 1];
        colInfo[colOrder[j]].pos = j;
      }
      colOrder[startCol + i] = col;
      colInfo[col].pos = startCol + i;
    }

    sccList = ALLOC(SccList_t, 1);
    sccList->nFuncs = array_n(rowScc);
    sccList->nVars = array_n(colScc);
    sccList->startFunc = startRow;
    sccList->lastFunc = startRow + sccList->nFuncs - 1;
    sccList->startVar = startCol;
    sccList->lastVar = startCol + sccList->nVars - 1;
    sccList->next = NIL(SccList_t);
    startRow += sccList->nFuncs;
    startCol += sccList->nVars;

    if (sccTailList)
      sccTailList->next = sccList;
    else
      sccHeadList = sccList;
    sccTailList = sccList;

    array_free(rowScc);
    array_free(colScc);

    if (option->mlpDebug) {
      CheckMatrix(xy, NIL(SccList_t), nActiveRows, nActiveCols,
                  rowInfo, colInfo, rowOrder, colOrder,
                  0, nActiveRows - 1, 0, nActiveCols - 1, 1);
    }
  }

  FREE(stack);
  FREE(stackFlag);
  FREE(rowFlag);
  FREE(colFlag);

  if (option->mlpSortScc && clusteredFlag == 0) {
    sccArray = array_alloc(SccList_t *, 0);
    sccList = sccHeadList;
    while (sccList) {
      array_insert_last(SccList_t *, sccArray, sccList);
      sccList = sccList->next;
    }

    if (option->mlpSortSccMode == 0) {
      if (option->mlpSortScc == 1)
        array_sort(sccArray, SccSortListIncreasingWithArea);
      else
        array_sort(sccArray, SccSortListDecreasingWithArea);
    } else if (option->mlpSortSccMode == 1) {
      if (option->mlpSortScc == 1)
        array_sort(sccArray, SccSortListIncreasingWithVars);
      else
        array_sort(sccArray, SccSortListDecreasingWithVars);
    } else {
      if (option->mlpSortScc == 1)
        array_sort(sccArray, SccSortListIncreasingWithRatio);
      else
        array_sort(sccArray, SccSortListDecreasingWithRatio);
    }

    newRowOrder = ALLOC(int, nActiveRows);
    newColOrder = ALLOC(int, nActiveCols);

    nCurRows = 0;
    nCurCols = 0;
    sccHeadList = NIL(SccList_t);
    for (i = 0; i < array_n(sccArray); i++) {
      sccList = array_fetch(SccList_t *, sccArray, i);
      if (sccHeadList)
        sccTailList->next = sccList;
      else
        sccHeadList = sccList;
      sccTailList = sccList;
      sccList->next = NIL(SccList_t);

      for (j = sccList->startFunc; j <= sccList->lastFunc; j++) {
        row = rowOrder[j];
        newRowOrder[nCurRows + j - sccList->startFunc] = rowOrder[j];
        rowInfo[row].pos = nCurRows + j - sccList->startFunc;
      }
      sccList->startFunc = nCurRows;
      sccList->lastFunc = nCurRows + sccList->nFuncs - 1;
      nCurRows += sccList->nFuncs;

      for (j = sccList->startVar; j <= sccList->lastVar; j++) {
        col = colOrder[j];
        newColOrder[nCurCols + j - sccList->startVar] = colOrder[j];
        colInfo[col].pos = nCurCols + j - sccList->startVar;
      }
      sccList->startVar = nCurCols;
      sccList->lastVar = nCurCols + sccList->nVars - 1;
      nCurCols += sccList->nVars;
    }

    memcpy(rowOrder, newRowOrder, sizeof(int) * nActiveRows);
    memcpy(colOrder, newColOrder, sizeof(int) * nActiveCols);

    if (option->mlpDebug) {
      CheckMatrix(xy, NIL(SccList_t), nActiveRows, nActiveCols,
                  rowInfo, colInfo, rowOrder, colOrder,
                  0, nActiveRows - 1, 0, nActiveCols - 1, 1);
    }

    FREE(newRowOrder);
    FREE(newColOrder);
    array_free(sccArray);
  }

  if (option->mlpVerbosity >= 2) {
    finalTime = util_cpu_time();
    fprintf(vis_stdout, "time for scc decomposition = %10g (#scc=%d)\n",
            (double)(finalTime - initialTime) / 1000.0, nGroups);
  }

  return(sccHeadList);
}


static int
SccSortListIncreasingWithVars(const void *p1, const void *p2)
{
  SccList_t **scc1 = (SccList_t **) p1;
  SccList_t **scc2 = (SccList_t **) p2;
  if ((*scc1)->nVars > (*scc2)->nVars)
    return(1);
  else if ((*scc1)->nVars < (*scc2)->nVars)
    return(-1);
  else
    return(0);
}


static int
SccSortListDecreasingWithVars(const void *p1, const void *p2)
{
  SccList_t **scc1 = (SccList_t **) p1;
  SccList_t **scc2 = (SccList_t **) p2;
  if ((*scc1)->nVars < (*scc2)->nVars)
    return(1);
  else if ((*scc1)->nVars > (*scc2)->nVars)
    return(-1);
  else
    return(0);
}


static int
SccSortListIncreasingWithArea(const void *p1, const void *p2)
{
  SccList_t **scc1 = (SccList_t **) p1;
  SccList_t **scc2 = (SccList_t **) p2;
  int   area1, area2;

  area1 = (*scc1)->nFuncs * (*scc1)->nVars;
  area2 = (*scc2)->nFuncs * (*scc2)->nVars;
  if (area1 > area2)
    return(1);
  else if (area1 < area2)
    return(-1);
  else
    return(0);
}


static int
SccSortListDecreasingWithArea(const void *p1, const void *p2)
{
  SccList_t **scc1 = (SccList_t **) p1;
  SccList_t **scc2 = (SccList_t **) p2;
  int   area1, area2;

  area1 = (*scc1)->nFuncs * (*scc1)->nVars;
  area2 = (*scc2)->nFuncs * (*scc2)->nVars;
  if (area1 < area2)
    return(1);
  else if (area1 > area2)
    return(-1);
  else
    return(0);
}


static int
SccSortListIncreasingWithRatio(const void *p1, const void *p2)
{
  SccList_t **scc1 = (SccList_t **) p1;
  SccList_t **scc2 = (SccList_t **) p2;
  float ratio1, ratio2;

  ratio1 = (float)((*scc1)->nVars) / (float)((*scc1)->nFuncs);
  ratio2 = (float)((*scc2)->nVars) / (float)((*scc2)->nFuncs);
  if (ratio1 > ratio2)
    return(1);
  else if (ratio1 < ratio2)
    return(-1);
  else
    return(0);
}


static int
SccSortListDecreasingWithRatio(const void *p1, const void *p2)
{
  SccList_t **scc1 = (SccList_t **) p1;
  SccList_t **scc2 = (SccList_t **) p2;
  float ratio1, ratio2;

  ratio1 = (float)((*scc1)->nVars) / (float)((*scc1)->nFuncs);
  ratio2 = (float)((*scc2)->nVars) / (float)((*scc2)->nFuncs);
  if (ratio1 < ratio2)
    return(1);
  else if (ratio1 > ratio2)
    return(-1);
  else
    return(0);
}


static int
SccSortRc(const void * p1, const void * p2)
{
  int *n1 = (int *) p1;
  int *n2 = (int *) p2;
  if (*n1 > *n2)
    return(1);
  else if (*n1 < *n2)
    return(-1);
  else
    return(0);
}


static void
FreeSccList(SccList_t *sccHeadList)
{
  SccList_t     *sccList, *sccNextList;

  sccList = sccHeadList;
  while (sccList) {
    sccNextList = sccList->next;
    FREE(sccList);
    sccList = sccNextList;
  }
}


static int
NumOfSccs(SccList_t *sccHeadList)
{
  SccList_t     *sccList;
  int           num = 0;

  sccList = sccHeadList;
  while (sccList) {
    num++;
    sccList = sccList->next;
  }
  return(num);
}


static void
BlockLowerTriangle(char **xy, int nRows, int nCols,
                   int nActiveRows, int nActiveCols,
                   SccList_t *sccList, int *rowOrder, int *colOrder,
                   RcInfo_t *rowInfo, RcInfo_t *colInfo,
                   ImgTrmOption_t *option)
{
  long          initialTime, finalTime;
  int           startFunc, lastFunc, startVar, lastVar;
  int           nLeftRows, nLeftCols;
  int           saveStartFunc, saveLastFunc;
  int           saveStartVar, saveLastVar;

  if (sccList) {
    if (sccList->nFuncs == 1 || sccList->nVars == 1)
      return;
    startFunc = sccList->startFunc;
    lastFunc = sccList->lastFunc;
    startVar = sccList->startVar;
    lastVar = sccList->lastVar;
  } else {
    startFunc = 0;
    lastFunc = nActiveRows - 1;
    startVar = 0;
    lastVar = nActiveCols - 1;
  }

  if (option->mlpVerbosity >= 2) {
    saveStartFunc = startFunc;
    saveLastFunc = lastFunc;
    saveStartVar = startVar;
    saveLastVar = lastVar;
  } else {
    /* To remove compiler warnings */
    saveStartFunc = 0;
    saveLastFunc = 0;
    saveStartVar = 0;
    saveLastVar = 0;
  }

  if (option->mlpDebug) {
    CheckMatrix(xy, sccList, nActiveRows, nActiveCols, rowInfo, colInfo,
                rowOrder, colOrder, startFunc, lastFunc, startVar, lastVar, 1);
  }

  if (option->mlpCsFirst) {
    /*
     * Find singleton columns and move each column to the right of the active
     * window and move its row to the bottom of the active window.
     */
    if (option->mlpSkipCs == 0) {
      if (option->mlpVerbosity >= 2)
        initialTime = util_cpu_time();
      else /* to remove uninitialized variables warning */
        initialTime = 0;
      FindAndMoveSingletonCols(xy, sccList, nActiveRows, nActiveCols,
                               &startFunc, &lastFunc, &startVar, &lastVar,
                               rowInfo, colInfo, rowOrder, colOrder, option);
      if (option->mlpVerbosity >= 2) {
        finalTime = util_cpu_time();
        fprintf(vis_stdout,
                "time for col singleton = %10g (#row=%d, #col=%d)\n",
                (double)(finalTime - initialTime) / 1000.0,
                saveLastFunc - lastFunc, saveLastVar - lastVar);
      }
    }

    /*
     * Find singleton rows and move each row to the top of the active window
     * and move its column to the left of the active window.
     */
    if (option->mlpSkipRs == 0) {
      if (option->mlpVerbosity >= 2)
        initialTime = util_cpu_time();
      else /* to remove uninitialized variables warning */
        initialTime = 0;
      FindAndMoveSingletonRows(xy, sccList, nActiveRows, nActiveCols,
                               &startFunc, &lastFunc, &startVar, &lastVar,
                               rowInfo, colInfo, rowOrder, colOrder, option);
      if (option->mlpVerbosity >= 2) {
        finalTime = util_cpu_time();
        fprintf(vis_stdout,
                "time for row singleton = %10g (#row=%d, #col=%d)\n",
                (double)(finalTime - initialTime) / 1000.0,
                startFunc - saveStartFunc, startVar - saveStartVar);
      }
    }
  } else {
    /*
     * Find singleton rows and move each row to the top of the active window
     * and move its column to the left of the active window.
     */
    if (option->mlpSkipRs == 0) {
      if (option->mlpVerbosity >= 2)
        initialTime = util_cpu_time();
      else /* to remove uninitialized variables warning */
        initialTime = 0;
      FindAndMoveSingletonRows(xy, sccList, nActiveRows, nActiveCols,
                               &startFunc, &lastFunc, &startVar, &lastVar,
                               rowInfo, colInfo, rowOrder, colOrder, option);
      if (option->mlpVerbosity >= 2) {
        finalTime = util_cpu_time();
        fprintf(vis_stdout,
                "time for row singleton = %10g (#row=%d, #col=%d)\n",
                (double)(finalTime - initialTime) / 1000.0,
                startFunc - saveStartFunc, startVar - saveStartVar);
      }
    }

    /*
     * Find singleton columns and move each column to the right of the active
     * window and move its row to the bottom of the active window.
     */
    if (option->mlpSkipCs == 0) {
      if (option->mlpVerbosity >= 2)
        initialTime = util_cpu_time();
      else /* to remove uninitialized variables warning */
        initialTime = 0;
      FindAndMoveSingletonCols(xy, sccList, nActiveRows, nActiveCols,
                               &startFunc, &lastFunc, &startVar, &lastVar,
                               rowInfo, colInfo, rowOrder, colOrder, option);
      if (option->mlpVerbosity >= 2) {
        finalTime = util_cpu_time();
        fprintf(vis_stdout,
                "time for col singleton = %10g (#row=%d, #col=%d)\n",
                (double)(finalTime - initialTime) / 1000.0,
                saveLastFunc - lastFunc, saveLastVar - lastVar);
      }
    }
  }

  if (option->mlpVerbosity >= 2) {
    initialTime = util_cpu_time();
    nLeftRows = lastFunc - startFunc + 1;
    nLeftCols = lastVar - startVar + 1;
  } else {
    /* To remove compiler warnings */
    initialTime = 0;
    nLeftRows = 0;
    nLeftCols = 0;
  }

  if (option->mlpRowBased) {
    MoveBestRows(xy, sccList, nActiveRows, nActiveCols, rowOrder, colOrder,
                 rowInfo, colInfo,
                 startFunc, lastFunc, startVar, lastVar, option);
  } else {
    MoveBestCols(xy, sccList, nRows, nCols, nActiveRows, nActiveCols,
                 rowOrder, colOrder, rowInfo, colInfo,
                 startFunc, lastFunc, startVar, lastVar, option);
  }

  if (option->mlpVerbosity >= 4 && option->mlpPrintMatrix) {
    PrintMatrix(xy, nActiveRows, nActiveCols, rowOrder, colOrder,
                rowInfo, colInfo, startFunc, lastFunc, startVar, lastVar);
  }

  if (option->mlpVerbosity >= 2) {
    finalTime = util_cpu_time();
    if (option->mlpRowBased) {
      fprintf(vis_stdout, "time for best row = %10g (#row=%d, #col=%d)\n",
                (double)(finalTime - initialTime) / 1000.0,
                nLeftRows - lastFunc + startFunc - 1,
                nLeftCols - lastVar + startVar - 1);
    } else {
      fprintf(vis_stdout, "time for best col = %10g (#row=%d, #col=%d)\n",
                (double)(finalTime - initialTime) / 1000.0,
                nLeftRows - lastFunc + startFunc - 1,
                nLeftCols - lastVar + startVar - 1);
    }
  }
}


static void
MoveBestRows(char **xy, SccList_t *sccList, int nRows, int nCols,
                int *rowOrder, int *colOrder,
                RcInfo_t *rowInfo, RcInfo_t *colInfo,
                int startFunc, int lastFunc, int startVar, int lastVar,
                ImgTrmOption_t *option)
{
  int           x, y, s1, t1, s2, t2;
  int           min, bestX;
  int           row, col, lastCol, moveRow, bestRow;
  int           regionX, overlap, nOverlaps, nMore;

  regionX = 0;
  bestX = -1;
  while (startFunc < lastFunc && startVar < lastVar) {
    min = nCols;
    overlap = 0;
    for (x = startFunc; x <= lastFunc; x++) {
      row = rowOrder[x];
      assert(rowInfo[row].lNum);
      if (overlap) {
        if (rowInfo[row].lNum < min) {
          s1 = colInfo[rowInfo[row].gMin].pos;
          if (s1 < regionX)
            s1 = regionX;
          nOverlaps = 0;
          for (y = s1; y < startVar; y++) {
            col = colOrder[y];
            if (xy[row][col]) {
              nOverlaps++;
              break;
            }
          }

          if (nOverlaps) {
            bestX = x;
            min = rowInfo[row].lNum;
          }
        }
      } else {
        s1 = colInfo[rowInfo[row].gMin].pos;
        if (s1 < regionX)
          s1 = regionX;
        nOverlaps = 0;
        for (y = s1; y < startVar; y++) {
          col = colOrder[y];
          if (xy[row][col]) {
            nOverlaps++;
            break;
          }
        }

        if (nOverlaps) {
          bestX = x;
          min = rowInfo[row].lNum;
          overlap = 1;
        } else {
          if (rowInfo[row].lNum < min) {
            bestX = x;
            min = rowInfo[row].lNum;
          }
        }
      }
    }

    assert(bestX != -1);
    if (!overlap)
      regionX = startFunc;

    row = rowOrder[bestX];
    while (1) {
      bestRow = row;
      s1 = colInfo[rowInfo[row].lMin].pos;
      t1 = colInfo[rowInfo[row].lMax].pos;
      min = rowInfo[row].lNum;
      for (x = startFunc; x <= lastFunc; x++) {
        moveRow = rowOrder[x];
        if (moveRow == row)
          continue;
        s2 = colInfo[rowInfo[moveRow].lMin].pos;
        t2 = colInfo[rowInfo[moveRow].lMax].pos;
        if (s2 < s1 || t2 > t1)
          continue;
        nMore = 0;
        nOverlaps = 0;
        for (y = s2; y <= t2; y++) {
          col = colOrder[y];
          if (xy[row][col]) {
            if (xy[moveRow][col])
              nOverlaps++;
          } else {
            if (xy[moveRow][col]) {
              nMore++;
              break;
            }
          }
        }
        if (nMore || nOverlaps == 0 || nOverlaps >= min)
          continue;
        min = nOverlaps;
        bestRow = moveRow;
      }

      lastCol = -1;
      s1 = colInfo[rowInfo[bestRow].lMin].pos;
      t1 = colInfo[rowInfo[bestRow].lMax].pos;
      for (y = s1; y <= t1; y++) {
        col = colOrder[y];
        if (xy[bestRow][col]) {
          if (option->mlpVerbosity >= 3) {
            fprintf(vis_stdout, "[%d] Moving Col %d(%d) to %d (best row)\n",
                    nMoves, y, col, startVar);
            fprintf(vis_stdout, "active window [%d-%d, %d-%d]\n",
                    startFunc, lastFunc, startVar, lastVar);
          }
          MoveColToLeft(xy, y, startFunc, lastFunc, startVar, lastVar,
                          rowInfo, colInfo, rowOrder, colOrder,
                          NIL(RcListInfo_t), option->mlpMethod);
          startVar++;
          lastCol = col;

          if (option->mlpVerbosity >= 4 && option->mlpPrintMatrix) {
            PrintMatrix(xy, nRows, nCols, rowOrder, colOrder, rowInfo, colInfo,
                        startFunc, lastFunc, startVar, lastVar);
          }
          if (option->mlpDebug) {
            CheckMatrix(xy, sccList, nRows, nCols, rowInfo, colInfo,
                        rowOrder, colOrder,
                        startFunc, lastFunc, startVar, lastVar, 1);
          }
        }
      }

      if (lastCol != -1) {
        s2 = rowInfo[colInfo[lastCol].lMin].pos;
        t2 = rowInfo[colInfo[lastCol].lMax].pos;
        for (x = s2; x <= t2; x++) {
          moveRow = rowOrder[x];
          if (xy[moveRow][lastCol] && rowInfo[moveRow].lNum == 0) {
            if (option->mlpVerbosity >= 3) {
              fprintf(vis_stdout, "[%d] Moving row %d(%d) to %d (best row)\n",
                        nMoves, x, moveRow, startFunc);
              fprintf(vis_stdout, "active window [%d-%d, %d-%d]\n",
                        startFunc, lastFunc, startVar, lastVar);
            }
            MoveRowToTop(xy, x, startFunc, lastFunc, startVar, lastVar,
                         rowInfo, colInfo, rowOrder, colOrder,
                         NIL(RcListInfo_t), option->mlpMethod);
            startFunc++;
            if (option->mlpVerbosity >= 4 && option->mlpPrintMatrix) {
              PrintMatrix(xy, nRows, nCols, rowOrder, colOrder,
                          rowInfo, colInfo,
                          startFunc, lastFunc, startVar, lastVar);
            }
            if (option->mlpDebug) {
              CheckMatrix(xy, sccList, nRows, nCols, rowInfo, colInfo,
                          rowOrder, colOrder,
                          startFunc, lastFunc, startVar, lastVar, 1);
            }
          }
        }
      }

      if (row == bestRow)
        break;
    }
  }
}


static void
MoveBestCols(char **xy, SccList_t *sccList, int nRows, int nCols,
                int nActiveRows, int nActiveCols,
                int *rowOrder, int *colOrder,
                RcInfo_t *rowInfo, RcInfo_t *colInfo,
                int startFunc, int lastFunc, int startVar, int lastVar,
                ImgTrmOption_t *option)
{
  int           x, y, i, s1, t1, s2, t2;
  char          *rowFlag, *colFlag;
  int           min, bestY, bestCol;
  int           row, col, tie, newTies, intersect = 0, maxIntersect, maxLength;
  int           minX, maxX;
  RcListInfo_t  *sortedRows, *sortedCols;
  RcList_t      *cur, *list;

  if (option->mlpSortedRows) {
    sortedRows = SortedListAlloc();
    for (x = startFunc; x <= lastFunc; x++) {
      row = rowOrder[x];
      SortedListInsert(sortedRows, row, row, rowInfo, 3);
    }

    if (option->mlpSortedCols)
      sortedCols = SortedListAlloc();
    else
      sortedCols = NIL(RcListInfo_t);

    if (option->mlpDebug)
      CheckSortedList(sortedRows, rowInfo, 3);
  } else {
    sortedRows = NIL(RcListInfo_t);
    sortedCols = NIL(RcListInfo_t);
  }

  rowFlag = ALLOC(char, nRows);
  if (sortedCols)
    memset(rowFlag, 0, sizeof(char) * nRows);
  colFlag = ALLOC(char, nCols);
  while (startFunc != lastFunc && startVar != lastVar) {
    min = nActiveCols;
    if (!sortedCols)
      memset(rowFlag, 0, sizeof(char) * nRows);
    tie = 0;
    minX = lastFunc;
    maxX = startFunc;
    if (sortedCols) {
      bestY = nActiveCols;
      bestCol = nActiveCols;
      if (sortedCols->num == 0) {
        if (sortedRows) {
          cur = sortedRows->head;
          while (cur) {
            row = cur->index;
            x = rowInfo[row].pos;
            if (rowInfo[row].lNum == 0) {
              cur = cur->next;
              continue;
            }
            if (rowInfo[row].lNum > min)
              break;
            rowFlag[row] = 1;
            s1 = colInfo[rowInfo[row].lMin].pos;
            t1 = colInfo[rowInfo[row].lMax].pos;
            for (y = s1; y <= t1; y++) {
              col = colOrder[y];
              if (!xy[row][col])
                continue;
              if (st_lookup(sortedCols->table, (char *)(long)col, &list)) {
                list->otherIndex++;
                SortedListMove(sortedCols, colInfo, col, 4);
              } else
                SortedListInsert(sortedCols, col, 1, colInfo, 4);
            }
            min = rowInfo[row].lNum;
            tie++;
            cur = cur->next;
          }
        }
      } else {
        if (sortedRows) {
          cur = sortedRows->head;
          while (cur) {
            row = cur->index;
            x = rowInfo[row].pos;
            if (rowInfo[row].lNum == 0) {
              cur = cur->next;
              continue;
            }
            if (rowInfo[row].lNum > min)
              break;
            min = rowInfo[row].lNum;
            tie++;
            cur = cur->next;
          }
        }
      }

      if (option->mlpDebug) {
        if (sortedRows)
          CheckSortedList(sortedRows, rowInfo, 3);
        CheckSortedList(sortedCols, colInfo, 4);
      }

      bestCol = sortedCols->head->index;
      bestY = colInfo[bestCol].pos;
      assert(bestY >= startVar && bestY <= lastVar);
      maxLength = colInfo[bestCol].lNum;
      intersect = sortedCols->head->otherIndex;
    } else {
      if (sortedRows) {
        if (option->mlpDebug)
          CheckSortedList(sortedRows, rowInfo, 3);

        cur = sortedRows->head;
        while (cur) {
          row = cur->index;
          x = rowInfo[row].pos;
          if (rowInfo[row].lNum == 0) {
            cur = cur->next;
            continue;
          }
          if (rowInfo[row].lNum < min) {
            min = rowInfo[row].lNum;
            tie = 1;
            memset(rowFlag, 0, sizeof(char) * nRows);
            rowFlag[x] = 1;
            minX = x;
            maxX = x;
          } else if (rowInfo[row].lNum == min) {
            tie++;
            rowFlag[x] = 1;
            if (x < minX)
              minX = x;
            if (x > maxX)
              maxX = x;
          } else
            break;
          cur = cur->next;
        }
      } else {
        for (x = startFunc; x <= lastFunc; x++) {
          row = rowOrder[x];
          if (rowInfo[row].lNum == 0)
            continue;
          if (rowInfo[row].lNum < min) {
            min = rowInfo[row].lNum;
            tie = 1;
            memset(rowFlag, 0, sizeof(char) * nRows);
            rowFlag[x] = 1;
            minX = x;
            maxX = x;
          } else if (rowInfo[row].lNum == min) {
            tie++;
            rowFlag[x] = 1;
            maxX = x;
          }
        }
      }

      if (tie == 0)
        break;

      maxLength = 0;
      maxIntersect = 0;
      memset(colFlag, 0, sizeof(char) * nCols);
      bestY = nActiveCols;
      bestCol = nActiveCols;
      for (x = minX; x <= maxX; x++) {
        if (!rowFlag[x])
          continue;
        row = rowOrder[x];
        s1 = colInfo[rowInfo[row].lMin].pos;
        t1 = colInfo[rowInfo[row].lMax].pos;
        for (y = s1; y <= t1; y++) {
          col = colOrder[y];
          if (xy[row][col]) {
            if (colFlag[y])
              continue;
            colFlag[y] = 1;
            intersect = 0;
            s2 = rowInfo[colInfo[col].lMin].pos;
            t2 = rowInfo[colInfo[col].lMax].pos;
            if (s2 < minX)
              s2 = minX;
            if (t2 > maxX)
              t2 = maxX;
            for (i = s2; i <= t2; i++) {
              if (xy[rowOrder[i]][col]) {
                if (rowFlag[i])
                  intersect++;
              }
            }
            if (intersect > maxIntersect ||
                (intersect == maxIntersect && colInfo[col].lNum > maxLength) ||
                (intersect == maxIntersect && colInfo[col].lNum == maxLength &&
                  y < bestY)) {
              bestY = y;
              bestCol = colOrder[bestY];
              maxLength = colInfo[bestCol].lNum;
              maxIntersect = intersect;
            }
          }
        }
      }
    }

    if (option->mlpVerbosity >= 3) {
      fprintf(vis_stdout,
"[%d] Moving Col %d(%d) to %d (best col, min=%d, intersect=%d, length=%d)\n",
             nMoves, bestY, bestCol, startVar, min, intersect, maxLength);
      fprintf(vis_stdout, "active window [%d-%d, %d-%d]\n",
             startFunc, lastFunc, startVar, lastVar);
    }
    MoveColToLeft(xy, bestY, startFunc, lastFunc, startVar, lastVar,
                  rowInfo, colInfo, rowOrder, colOrder,
                  NIL(RcListInfo_t), option->mlpMethod);
    startVar++;
    if (option->mlpVerbosity >= 4 && option->mlpPrintMatrix) {
      PrintMatrix(xy, nActiveRows, nActiveCols, rowOrder, colOrder,
                  rowInfo, colInfo, startFunc, lastFunc, startVar, lastVar);
    }
    if (option->mlpDebug) {
      CheckMatrix(xy, sccList, nActiveRows, nActiveCols, rowInfo, colInfo,
                  rowOrder, colOrder,
                  startFunc, lastFunc, startVar, lastVar, 1);
    }

    if (sortedCols)
      SortedListDelete(sortedCols, bestCol);

    if (sortedRows) {
      if (min == 1) {
        cur = sortedRows->head;
        while (cur) {
          row = cur->index;
          if (rowInfo[row].lNum > min)
            break;
          else if (rowInfo[row].lNum == min) {
            cur = cur->next;
            continue;
          }
          x = rowInfo[row].pos;
          if (option->mlpVerbosity >= 3) {
            fprintf(vis_stdout, "[%d] Moving row %d(%d) to %d (empty row)\n",
                   nMoves, x, row, startFunc);
            fprintf(vis_stdout, "active window [%d-%d, %d-%d]\n",
                   startFunc, lastFunc, startVar, lastVar);
          }
          MoveRowToTop(xy, x, startFunc, lastFunc, startVar, lastVar,
                       rowInfo, colInfo, rowOrder, colOrder,
                       NIL(RcListInfo_t), option->mlpMethod);
          startFunc++;
          if (option->mlpVerbosity >= 4 && option->mlpPrintMatrix) {
            PrintMatrix(xy, nActiveRows, nActiveCols, rowOrder, colOrder,
                        rowInfo, colInfo,
                        startFunc, lastFunc, startVar, lastVar);
          }
          if (option->mlpDebug) {
            CheckMatrix(xy, sccList, nActiveRows, nActiveCols,
                        rowInfo, colInfo, rowOrder, colOrder,
                        startFunc, lastFunc, startVar, lastVar, 1);
          }
          cur = cur->next;
          SortedListDelete(sortedRows, row);
        }
      }

      /* Update sortedRows */
      /*
      cur = sortedRows->head;
      while (cur) {
        row = cur->index;
        if (xy[row][bestCol]) {
          assert(rowInfo[row].lNum > 0);
          SortedListMove(sortedRows, rowInfo, row, 3);
        }
        cur = cur->next;
      }
      if (option->mlpDebug)
        CheckSortedList(sortedRows, rowInfo, 3);
      */
      s1 = rowInfo[colInfo[bestCol].lMin].pos;
      t1 = rowInfo[colInfo[bestCol].lMax].pos;
      cur = sortedRows->head;
      while (cur) {
        row = cur->index;
        x = rowInfo[row].pos;
        if (x < s1) {
          cur = cur->next;
          continue;
        }
        if (xy[row][bestCol] && rowInfo[row].lNum > 0) {
          SortedListMove(sortedRows, rowInfo, row, 3);
          if (sortedCols && min == 1) {
            if (rowInfo[row].lNum == 1 && rowFlag[row] == 0) {
              rowFlag[row] = 1;
              s2 = colInfo[rowInfo[row].lMin].pos;
              t2 = colInfo[rowInfo[row].lMax].pos;
              for (y = s2; y <= t2; y++) {
                col = colOrder[y];
                if (!xy[row][col])
                  continue;
                if (st_lookup(sortedCols->table, (char *)(long)col, &list)) {
                  list->otherIndex++;
                  SortedListMove(sortedCols, colInfo, col, 4);
                } else
                  SortedListInsert(sortedCols, col, 1, colInfo, 4);
              }
            }
          }
        }
        cur = cur->next;
      }

      newTies = 0;
      if (sortedCols && min > 1) {
        min = nActiveCols;
        cur = sortedRows->head;
        while (cur) {
          row = cur->index;
          x = rowInfo[row].pos;
          if (rowInfo[row].lNum == 0) {
            cur = cur->next;
            continue;
          }
          if (rowInfo[row].lNum > min)
            break;
          min = rowInfo[row].lNum;
          newTies++;
          cur = cur->next;
        }

        if (newTies != tie) {
          SortedListFree(sortedCols);
          sortedCols = SortedListAlloc();

          cur = sortedRows->head;
          while (cur) {
            row = cur->index;
            x = rowInfo[row].pos;
            if (rowInfo[row].lNum == 0) {
              cur = cur->next;
              continue;
            }
            if (rowInfo[row].lNum > min)
              break;
            rowFlag[row] = 1;
            s1 = colInfo[rowInfo[row].lMin].pos;
            t1 = colInfo[rowInfo[row].lMax].pos;
            for (y = s1; y <= t1; y++) {
              col = colOrder[y];
              if (!xy[row][col])
                continue;
              if (st_lookup(sortedCols->table, (char *)(long)col, &list)) {
                list->otherIndex++;
                SortedListMove(sortedCols, colInfo, col, 4);
              } else
                SortedListInsert(sortedCols, col, 1, colInfo, 4);
            }
            cur = cur->next;
          }
        }
      }
    } else {
      if (min == 1) {
        s1 = rowInfo[colInfo[bestCol].lMin].pos;
        t1 = rowInfo[colInfo[bestCol].lMax].pos;
        for (x = s1; x <= t1; x++) {
          row = rowOrder[x];
          if (xy[row][bestCol] && rowInfo[row].lNum == 0) {
            if (option->mlpVerbosity >= 3) {
              fprintf(vis_stdout, "[%d] Moving row %d(%d) to %d (empty row)\n",
                     nMoves, x, row, startFunc);
              fprintf(vis_stdout, "active window [%d-%d, %d-%d]\n",
                     startFunc, lastFunc, startVar, lastVar);
            }
            MoveRowToTop(xy, x, startFunc, lastFunc, startVar, lastVar,
                         rowInfo, colInfo, rowOrder, colOrder,
                         NIL(RcListInfo_t), option->mlpMethod);
            startFunc++;
            if (option->mlpVerbosity >= 4 && option->mlpPrintMatrix) {
              PrintMatrix(xy, nActiveRows, nActiveCols, rowOrder, colOrder,
                          rowInfo, colInfo,
                          startFunc, lastFunc, startVar, lastVar);
            }
            if (option->mlpDebug) {
              CheckMatrix(xy, sccList, nActiveRows, nActiveCols,
                          rowInfo, colInfo, rowOrder, colOrder,
                          startFunc, lastFunc, startVar, lastVar, 1);
            }
          }
        }
      }
    }
  }

  if (sortedRows) {
    cur = sortedRows->head;
    while (cur) {
      row = cur->index;
      x = rowInfo[row].pos;
      if (option->mlpVerbosity >= 3) {
        fprintf(vis_stdout, "[%d] Moving row %d(%d) to %d (bandwidth)\n",
               nMoves, x, row, startFunc);
        fprintf(vis_stdout, "active window [%d-%d, %d-%d]\n",
               startFunc, lastFunc, startVar, lastVar);
      }
      MoveRowToTop(xy, x, startFunc, lastFunc, startVar, lastVar,
                   rowInfo, colInfo, rowOrder, colOrder,
                   NIL(RcListInfo_t), option->mlpMethod);
      startFunc++;
      if (option->mlpVerbosity >= 4 && option->mlpPrintMatrix) {
        PrintMatrix(xy, nActiveRows, nActiveCols, rowOrder, colOrder,
                    rowInfo, colInfo,
                    startFunc, lastFunc, startVar, lastVar);
      }
      if (option->mlpDebug) {
        CheckMatrix(xy, sccList, nActiveRows, nActiveCols, rowInfo, colInfo,
                    rowOrder, colOrder,
                    startFunc, lastFunc, startVar, lastVar, 1);
      }
      cur = cur->next;
      SortedListDelete(sortedRows, row);
    }
    SortedListFree(sortedRows);
  }

  if (sortedCols) {
    cur = sortedCols->head;
    while (cur) {
      col = cur->index;
      cur = cur->next;
      SortedListDelete(sortedCols, col);
    }
    SortedListFree(sortedCols);
  }

  FREE(rowFlag);
  FREE(colFlag);
}


static void
MlpPostProcess(char **xy, SccList_t *sccList, int nVars, int nRows, int nCols,
                RcInfo_t *rowInfo, RcInfo_t *colInfo,
                int *rowOrder, int *colOrder,
                Img_DirectionType direction, ImgTrmOption_t *option)
{
  int           x, y, i, j, r, c, row, col, otherRow, otherCol, tmpRow, tmpCol;
  int           nqv1, nqv2, niv1, niv2, benefit1, benefit2;
  int           s1, t1, s2, t2, s3, t3, s4, t4;
  long          initialTime, finalTime;
  int           nSwaps;
  float         lambda1, lambda2, lambda3;
  int           startFunc, lastFunc;

  if (sccList) {
    startFunc = sccList->startFunc;
    lastFunc = sccList->lastFunc;
  } else {
    startFunc = 0;
    lastFunc = nRows - 1;
  }

  if (option->mlpVerbosity)
    initialTime = util_cpu_time();
  else
    initialTime = 0;

  nSwaps = 0;
  for (x = lastFunc - 1; x >= startFunc; x--) {
    row = rowOrder[x];
    if (rowInfo[row].data.row.nsVarBddArray)
      niv1 = array_n(rowInfo[row].data.row.nsVarBddArray);
    else
      niv1 = 0;

    col = rowInfo[row].gMin;
    if (rowInfo[row].gNum == 1 &&
        colInfo[col].gNum > 1 && colInfo[col].gMin == x) {
      t1 = rowInfo[colInfo[col].gMax].pos;
      for (i = x + 1; i <= t1; i++) {
        otherRow = rowOrder[i];
        if (xy[otherRow][col]) {
          if (x < i - 1) {
            for (j = x; j < i - 1; j++) {
              rowOrder[j] = rowOrder[j + 1];
              rowInfo[rowOrder[j]].pos = j;
            }
            rowOrder[i - 1] = row;
            rowInfo[row].pos = i - 1;

            if (rowInfo[otherRow].gNum == 1)
              break;

            s2 = colInfo[rowInfo[otherRow].gMin].pos;
            t2 = colInfo[rowInfo[otherRow].gMax].pos;
            for (j = t2; j >= s2; j--) {
              otherCol = colOrder[j];
              if (colInfo[otherCol].gMin == otherRow)
                t2 = j - 1;
              else
                break;
            }

            s3 = rowInfo[colInfo[col].gMin].pos;
            t3 = rowInfo[colInfo[col].gMax].pos;
            for (r = s3; r <= t3; r++) {
              tmpRow = rowOrder[r];
              if (col == rowInfo[tmpRow].gMin) {
                s4 = colInfo[rowInfo[tmpRow].gMin].pos;
                t4 = colInfo[rowInfo[tmpRow].gMax].pos;
                if (t4 > t2)
                  t4 = t2;
                for (c = s4 + 1; c <= t4; c++) {
                  tmpCol = colOrder[c];
                  if (xy[tmpRow][tmpCol]) {
                    rowInfo[tmpRow].gMin = tmpCol;
                    rowInfo[tmpRow].lMin = tmpCol;
                    break;
                  }
                }
              }

              if (xy[tmpRow][col] && t2 >= colInfo[rowInfo[tmpRow].gMax].pos) {
                rowInfo[tmpRow].gMax = col;
                rowInfo[tmpRow].lMax = col;
              }
            }

            y = colInfo[col].pos;
            if (y < t2) {
              for (j = y; j < t2; j++) {
                colOrder[j] = colOrder[j + 1];
                colInfo[colOrder[j]].pos = j;
              }
              colOrder[t2] = col;
              colInfo[col].pos = t2;
            }

            if (option->mlpVerbosity >= 4 && option->mlpPrintMatrix) {
              PrintMatrix(xy, nRows, nCols, rowOrder, colOrder,
                          rowInfo, colInfo, 0, nRows - 1, 0, nCols - 1);
            }
            if (option->mlpDebug) {
              CheckMatrix(xy, sccList, nRows, nCols, rowInfo, colInfo,
                          rowOrder, colOrder, 0, nRows - 1, 0, nCols - 1, 1);
            }

            nSwaps++;
          }
          break;
        }
      }
      continue;
    }

    for (i = x + 1; i <= lastFunc; i++) {
      otherRow = rowOrder[i];
      if (rowInfo[otherRow].data.row.nsVarBddArray)
        niv2 = array_n(rowInfo[otherRow].data.row.nsVarBddArray);
      else
        niv2 = 0;

      nqv1 = 0;
      s1 = colInfo[rowInfo[row].gMin].pos;
      t1 = colInfo[rowInfo[row].gMax].pos;
      for (y = t1; y >= s1; y--) {
        col = colOrder[y];
        if (colInfo[col].gMin == row) {
          if (!xy[otherRow][col])
            nqv1++;
        } else
          break;
      }
      benefit1 = nqv1 - niv1;

      nqv2 = 0;
      s2 = colInfo[rowInfo[otherRow].gMin].pos;
      t2 = colInfo[rowInfo[otherRow].gMax].pos;
      for (y = t2; y >= s2; y--) {
        col = colOrder[y];
        if (colInfo[col].gMin == otherRow)
          nqv2++;
        else
          break;
      }
      benefit2 = nqv2 - niv2;

      if (benefit1 > benefit2 ||
          (benefit1 == benefit2 &&
           (rowInfo[row].gNum > rowInfo[otherRow].gNum ||
            (rowInfo[row].gNum == rowInfo[otherRow].gNum &&
             bdd_size(rowInfo[row].data.row.func) <
             bdd_size(rowInfo[otherRow].data.row.func))))) {
        nqv1 = 0;
        for (y = t1; y >= s1; y--) {
          if (t1 > t2)
            break;
          col = colOrder[y];
          if (colInfo[col].gMin == row) {
            if (!xy[otherRow][col]) {
              s3 = rowInfo[colInfo[col].gMin].pos;
              t3 = rowInfo[colInfo[col].gMax].pos;
              for (r = s3; r <= t3; r++) {
                tmpRow = rowOrder[r];
                if (col == rowInfo[tmpRow].gMin) {
                  s4 = colInfo[rowInfo[tmpRow].gMin].pos;
                  t4 = colInfo[rowInfo[tmpRow].gMax].pos;
                  if (t4 > t2 - nqv1)
                    t4 = t2 - nqv1;
                  for (c = s4 + 1; c <= t4; c++) {
                    tmpCol = colOrder[c];
                    if (xy[tmpRow][tmpCol] && colInfo[tmpCol].gMin != row) {
                      rowInfo[tmpRow].gMin = tmpCol;
                      rowInfo[tmpRow].lMin = tmpCol;
                      break;
                    }
                  }
                }

                if (nqv1 == 0 && xy[tmpRow][col] &&
                    t2 - nqv1 >= colInfo[rowInfo[tmpRow].gMax].pos) {
                  rowInfo[tmpRow].gMax = col;
                  rowInfo[tmpRow].lMax = col;
                }
              }

              for (j = y; j < t2 - nqv1; j++) {
                colOrder[j] = colOrder[j + 1];
                colInfo[colOrder[j]].pos = j;
              }
              colOrder[t2 - nqv1] = col;
              colInfo[col].pos = t2 - nqv1;
              nqv1++;
            } else {
              colInfo[col].gMin = otherRow;
              colInfo[col].lMin = otherRow;
              if (colInfo[col].gMax == otherRow) {
                colInfo[col].gMax = row;
                colInfo[col].lMax = row;
              }
            }
          } else {
            if (colInfo[col].gMax == otherRow && xy[row][col]) {
              colInfo[col].gMax = row;
              colInfo[col].lMax = row;
            }
          }
        }

        rowOrder[i - 1] = otherRow;
        rowInfo[otherRow].pos = i - 1;
        rowOrder[i] = row;
        rowInfo[row].pos = i;
        nSwaps++;

        if (option->mlpVerbosity >= 4 && option->mlpPrintMatrix) {
          PrintMatrix(xy, nRows, nCols, rowOrder, colOrder,
                        rowInfo, colInfo, 0, nRows - 1, 0, nCols - 1);
        }
        if (option->mlpDebug) {
          CheckMatrix(xy, sccList, nRows, nCols, rowInfo, colInfo,
                        rowOrder, colOrder, 0, nRows - 1, 0, nCols - 1, 1);
        }
      } else
        break;
    }
  }

  if (option->mlpVerbosity) {
    if (option->mlpVerbosity >= 2) {
      lambda1 = ComputeLambdaMlp(xy, nVars, nRows, nCols, rowInfo, colInfo,
                                 rowOrder, colOrder, direction, 0, 0, option);
      lambda2 = ComputeLambdaMlp(xy, nVars, nRows, nCols, rowInfo, colInfo,
                                 rowOrder, colOrder, direction, 1, 0, option);
      lambda3 = ComputeLambdaMlp(xy, nVars, nRows, nCols, rowInfo, colInfo,
                                 rowOrder, colOrder, direction, 2, 0, option);
      fprintf(vis_stdout, "Optimized Lambda = %f (%f, %f)\n",
              lambda1, lambda2, lambda3);
    }
    finalTime = util_cpu_time();
    fprintf(vis_stdout, "time for postprocessing = %10g (#swap=%d)\n",
           (double)(finalTime - initialTime) / 1000.0, nSwaps);
  }
}


static float
ComputeLambdaMlp(char **xy, int nVars, int nRows, int nCols,
                RcInfo_t *rowInfo, RcInfo_t *colInfo,
                int *rowOrder, int *colOrder,
                Img_DirectionType direction,
                int mode, int asIs, ImgTrmOption_t *option)
{
  int           lifetime, area;
  int           highest, lowest;
  int           x, y, row, col;
  int           nIntroducedVars, nYVars;
  float         lambda;

  lifetime = 0;

  if (direction == Img_Forward_c) {
    for (y = 0; y < nCols; y++) {
      col = colOrder[y];
      lowest = rowInfo[colInfo[col].gMin].pos;
      if (mode == 1) { /* active lifetime (alpha) */
        highest = rowInfo[colInfo[col].gMax].pos;
        lifetime += highest - lowest + 1;
      } else { /* total lifetime (lambda) */
        lifetime += nRows - lowest;
      }
    }
  } else {
    for (y = 0; y < nCols; y++) {
      col = colOrder[y];
      lowest = rowInfo[colInfo[col].gMin].pos;
      if (mode == 1) { /* active lifetime (alpha) */
        highest = rowInfo[colInfo[col].gMax].pos;
        lifetime += highest - lowest + 1;
      } else { /* total lifetime (lambda) */
        if (colInfo[col].data.col.type == 2) /* primary input variables */
          lifetime += lowest + 1;
        else { /* present state variables */
          if (mode == 2)
            lifetime += nRows - lowest;
        }
      }
    }
  }

  /* total lifetime (lambda) with introducedVars */
  nIntroducedVars = 0;;
  if (mode == 2 || (mode == 0 && direction == Img_Backward_c)) {
    if (asIs) { /* contains next state variables */
      for (x = 0; x < nRows; x++) {
        row = rowOrder[x];
        nYVars = array_n(rowInfo[row].data.row.nsVarBddArray);
        nIntroducedVars += nYVars;
        lifetime += (x + 1) * nYVars;
      }
    } else { /* does not contain next state variables */
      nIntroducedVars = nRows;
      lifetime += nIntroducedVars * (nIntroducedVars + 1) / 2;
    }
  }

  if (option->mlpLambdaMode == 0)
    area = nRows * (nCols + nIntroducedVars);
  else
    area = nRows * (nVars + nIntroducedVars);
  lambda = (float)lifetime / (float)area;
  return(lambda);
}


static void
FindAndMoveSingletonRows(char **xy, SccList_t *sccList, int nRows, int nCols,
        int *startFunc, int *lastFunc, int *startVar, int *lastVar,
        RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder,
        ImgTrmOption_t *option)
{
  int           x, row, col, add, found;
  RcListInfo_t  *candList;

  if (option->mlpMethod == 0) {
    while (*startFunc != *lastFunc && *startVar != *lastVar) {
      found = 0;
      for (x = *startFunc; x <= *lastFunc; x++) {
        row = rowOrder[x];
        if (rowInfo[row].lNum == 1) {
          found = 1;
          add = MoveSingletonRow(xy, sccList, nRows, nCols, x,
                                startFunc, lastFunc, startVar, lastVar,
                                rowInfo, colInfo, rowOrder, colOrder,
                                NIL(RcListInfo_t), option);
          if (*startFunc == *lastFunc || *startVar == *lastVar) {
            found = 0;
            break;
          }
          x += add;
        }
      }
      if (!found)
        break;
    }
  } else {
    if (*startFunc != *lastFunc && *startVar != *lastVar) {
      candList = SortedListAlloc();
      for (x = *startFunc; x <= *lastFunc; x++) {
        row = rowOrder[x];
        if (rowInfo[row].lNum == 1) {
          col = rowInfo[row].lMin;
          SortedListInsert(candList, row, col, colInfo, option->mlpMethod);
        }
      }
      if (option->mlpDebug)
        CheckSortedList(candList, colInfo, option->mlpMethod);
      while (candList->cur) {
        candList->curIndex = candList->cur->index;
        x = rowInfo[candList->curIndex].pos;
        MoveSingletonRow(xy, sccList, nRows, nCols, x,
                        startFunc, lastFunc, startVar, lastVar,
                        rowInfo, colInfo, rowOrder, colOrder,
                        candList, option);
        if (option->mlpDebug)
          CheckSortedList(candList, colInfo, option->mlpMethod);
        if (*startFunc == *lastFunc || *startVar == *lastVar) {
          break;
        }
      }
      if (option->mlpVerbosity) {
        fprintf(vis_stdout, "max number of elements in list = %d\n",
                candList->maxNum);
      }
      SortedListFree(candList);
    }
  }
}


static int
MoveSingletonRow(char **xy, SccList_t *sccList, int nRows, int nCols, int x,
        int *startFunc, int *lastFunc, int *startVar, int *lastVar,
        RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder,
        RcListInfo_t *candList, ImgTrmOption_t *option)
{
  int   y, i, add;
  int   row, col, otherRow;
  int   startRow, lastRow, startCol, lastCol;

  startRow = *startFunc;
  startCol = *startVar;
  lastRow = *lastFunc;
  lastCol = *lastVar;

  row = rowOrder[x];
  col = rowInfo[row].lMin;
  y = colInfo[col].pos;

  if (option->mlpVerbosity >= 3) {
    fprintf(vis_stdout, "[%d] Moving Col %d(%d) to %d (row singleton)\n",
          nMoves, y, col, startCol);
    fprintf(vis_stdout, "active window [%d-%d, %d-%d]\n",
           startRow, lastRow, startCol, lastCol);
  }
  MoveColToLeft(xy, y, startRow, lastRow, startCol, lastCol,
                rowInfo, colInfo, rowOrder, colOrder, candList,
                option->mlpMethod);
  startCol++;
  if (option->mlpVerbosity >= 4 && option->mlpPrintMatrix) {
    PrintMatrix(xy, nRows, nCols, rowOrder, colOrder, rowInfo, colInfo,
                startRow, lastRow, startCol, lastCol);
  }
  if (option->mlpDebug) {
    CheckMatrix(xy, sccList, nRows, nCols, rowInfo, colInfo, rowOrder, colOrder,
                startRow, lastRow, startCol, lastCol, 1);
  }

  add = 0;
  for (i = rowInfo[colInfo[col].lMin].pos;
       colInfo[col].lNum && i <= rowInfo[colInfo[col].lMax].pos; i++) {
    otherRow = rowOrder[i];
    if (!xy[otherRow][col])
      continue;
    if (rowInfo[otherRow].lNum == 0) {
      if (option->mlpVerbosity >= 3) {
        fprintf(vis_stdout, "[%d] Moving row %d(%d) to %d (row singleton)\n",
                nMoves, i, otherRow, startRow);
        fprintf(vis_stdout, "active window [%d-%d, %d-%d]\n",
               startRow, lastRow, startCol, lastCol);
      }
      MoveRowToTop(xy, i, startRow, lastRow, startCol, lastCol,
                   rowInfo, colInfo, rowOrder, colOrder,
                   NIL(RcListInfo_t), option->mlpMethod);
      startRow++;
      if (i > x)
        add++;
      if (option->mlpVerbosity >= 4 && option->mlpPrintMatrix) {
        PrintMatrix(xy, nRows, nCols, rowOrder, colOrder, rowInfo, colInfo,
                    startRow, lastRow, startCol, lastCol);
      }
      if (option->mlpDebug) {
        CheckMatrix(xy, sccList, nRows, nCols,
                    rowInfo, colInfo, rowOrder, colOrder,
                    startRow, lastRow, startCol, lastCol, 1);
      }
    }
  }

  *startFunc = startRow;
  *startVar = startCol;
  return(add);
}


static void
FindAndMoveSingletonCols(char **xy, SccList_t *sccList, int nRows, int nCols,
        int *startFunc, int *lastFunc, int *startVar, int *lastVar,
        RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder,
        ImgTrmOption_t *option)
{
  int           y, row, col, add, found;
  RcListInfo_t  *candList;

  if (option->mlpMethod == 0) {
    while (*startFunc != *lastFunc && *startVar != *lastVar) {
      found = 0;
      for (y = *startVar; y <= *lastVar; y++) {
        col = colOrder[y];
        if (colInfo[col].lNum == 1) {
          found = 1;
          add = MoveSingletonCol(xy, sccList, nRows, nCols, y,
                                startFunc, lastFunc, startVar, lastVar,
                                rowInfo, colInfo, rowOrder, colOrder,
                                NIL(RcListInfo_t), option);
          if (*startFunc == *lastFunc || *startVar == *lastVar) {
            found = 0;
            break;
          }
          y -= add;
        }
      }
      if (!found)
        break;
    }
  } else {
    if (*startFunc != *lastFunc && *startVar != *lastVar) {
      candList = SortedListAlloc();
      for (y = *startVar; y <= *lastVar; y++) {
        col = colOrder[y];
        if (colInfo[col].lNum == 1) {
          row = colInfo[col].lMin;
          SortedListInsert(candList, col, row, rowInfo, option->mlpMethod);
        }
      }
      if (option->mlpDebug)
        CheckSortedList(candList, rowInfo, option->mlpMethod);
      while (candList->cur) {
        candList->curIndex = candList->cur->index;
        y = colInfo[candList->curIndex].pos;
        MoveSingletonCol(xy, sccList, nRows, nCols, y,
                        startFunc, lastFunc, startVar, lastVar,
                        rowInfo, colInfo, rowOrder, colOrder, candList, option);
        if (option->mlpDebug)
          CheckSortedList(candList, rowInfo, option->mlpMethod);
        if (*startFunc == *lastFunc || *startVar == *lastVar) {
          break;
        }
      }
      if (option->mlpVerbosity >= 2) {
        fprintf(vis_stdout, "max number of elements in list = %d\n",
                candList->maxNum);
      }
      SortedListFree(candList);
    }
  }
}


static int
MoveSingletonCol(char **xy, SccList_t *sccList, int nRows, int nCols, int y,
        int *startFunc, int *lastFunc, int *startVar, int *lastVar,
        RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder,
        RcListInfo_t *candList, ImgTrmOption_t *option)
{
  int   x, j, add;
  int   row, col, otherCol;
  int   startRow, lastRow, startCol, lastCol;

  startRow = *startFunc;
  startCol = *startVar;
  lastRow = *lastFunc;
  lastCol = *lastVar;

  col = colOrder[y];
  row = colInfo[col].lMin;
  x = rowInfo[row].pos;

  if (option->mlpVerbosity >= 3) {
    fprintf(vis_stdout, "[%d] Moving Row %d(%d) to %d (col singleton)\n",
          nMoves, x, row, lastRow);
    fprintf(vis_stdout, "active window [%d-%d, %d-%d]\n",
           startRow, lastRow, startCol, lastCol);
  }
  MoveRowToBottom(xy, x, startRow, lastRow, startCol, lastCol,
                  rowInfo, colInfo, rowOrder, colOrder, candList,
                  option->mlpMethod);
  lastRow--;
  if (option->mlpVerbosity >= 4 && option->mlpPrintMatrix) {
    PrintMatrix(xy, nRows, nCols, rowOrder, colOrder, rowInfo, colInfo,
                startRow, lastRow, startCol, lastCol);
  }
  if (option->mlpDebug) {
    CheckMatrix(xy, sccList, nRows, nCols, rowInfo, colInfo, rowOrder, colOrder,
                startRow, lastRow, startCol, lastCol, 1);
  }

  add = 0;
  for (j = colInfo[rowInfo[row].lMin].pos;
       rowInfo[row].lNum && j <= colInfo[rowInfo[row].lMax].pos; j++) {
    otherCol = colOrder[j];
    if (!xy[row][otherCol])
      continue;
    if (colInfo[otherCol].lNum == 0) {
      if (option->mlpVerbosity >= 3) {
        fprintf(vis_stdout, "[%d] Moving Col %d(%d) to %d (col singleton)\n",
                nMoves, j, otherCol, lastCol);
        fprintf(vis_stdout, "active window [%d-%d, %d-%d]\n",
                startRow, lastRow, startCol, lastCol);
      }
      MoveColToRight(xy, j, startRow, lastRow, startCol, lastCol,
                     rowInfo, colInfo, rowOrder, colOrder,
                     NIL(RcListInfo_t), option->mlpMethod);
      lastCol--;
      j--;
      if (otherCol <= col)
        add++;
      if (option->mlpVerbosity >= 4 && option->mlpPrintMatrix) {
        PrintMatrix(xy, nRows, nCols, rowOrder, colOrder, rowInfo, colInfo,
                    startRow, lastRow, startCol, lastCol);
      }
      if (option->mlpDebug) {
        CheckMatrix(xy, sccList, nRows, nCols,
                    rowInfo, colInfo, rowOrder, colOrder,
                    startRow, lastRow, startCol, lastCol, 1);
      }
    }
  }

  *lastFunc = lastRow;
  *lastVar = lastCol;
  return(add);
}


static void
MoveRowToTop(char **xy, int x, int startFunc, int lastFunc,
             int startVar, int lastVar, RcInfo_t *rowInfo,
             RcInfo_t *colInfo, int *rowOrder, int *colOrder,
             RcListInfo_t *candList, int method)
{
  int   i, j, s, t;
  int   row, col, tmpRow;

  assert(x >= startFunc && x <= lastFunc);

  row = rowOrder[x];
  if (x != startFunc) {
    /* Change row order */
    for (i = x; i > startFunc; i--) {
      rowOrder[i] = rowOrder[i - 1];
      rowInfo[rowOrder[i]].pos = i;
    }
    rowOrder[startFunc] = row;
    rowInfo[row].pos = startFunc;
  }

  s = colInfo[rowInfo[row].gMin].pos;
  t = colInfo[rowInfo[row].gMax].pos;
  for (j = s; j <= t; j++) {
    col = colOrder[j];
    if (!xy[row][col])
      continue;
    colInfo[col].prevG = startFunc;
    if (colInfo[col].lNum == 1) {
      colInfo[col].lMin = row;
      colInfo[col].lMax = row;
    } else {
      if (row == colInfo[col].lMin) {
        for (i = x + 1; i <= rowInfo[colInfo[col].lMax].pos; i++) {
          tmpRow = rowOrder[i];
          if (xy[tmpRow][col]) {
            colInfo[col].lMin = tmpRow;
            break;
          }
        }
      } else if (row == colInfo[col].lMax) {
        for (i = x; i >= rowInfo[colInfo[col].lMin].pos; i--) {
          tmpRow = rowOrder[i];
          if (xy[tmpRow][col]) {
            colInfo[col].lMax = tmpRow;
            break;
          }
        }
      }
    }

    if (colInfo[col].gNum == 1) {
      colInfo[col].gMin = row;
      colInfo[col].gMax = row;
    } else {
      if (row == colInfo[col].gMin)
        colInfo[col].gMin = row;
      else if (row == colInfo[col].gMax) {
        colInfo[col].gMax = colInfo[col].lMax;
        if (startFunc <= rowInfo[colInfo[col].gMin].pos)
          colInfo[col].gMin = row;
      } else if (startFunc <= rowInfo[colInfo[col].gMin].pos)
        colInfo[col].gMin = row;
    }

    if (rowInfo[row].lNum) {
      if (j == rowInfo[row].prevG)
        j = colInfo[rowInfo[row].lMin].pos - 1;
      else if (j == colInfo[rowInfo[row].lMax].pos)
        j = rowInfo[row].nextG - 1;
    }
  }

  s = colInfo[rowInfo[row].lMin].pos;
  t = colInfo[rowInfo[row].lMax].pos;
  for (j = s; j <= t; j++) {
    col = colOrder[j];
    if (xy[row][col]) {
      colInfo[col].lNum--;
      if (candList) {
        if (colInfo[col].lNum == 0)
          SortedListDelete(candList, col);
        else if (colInfo[col].lNum == 1)
          SortedListInsert(candList, col, row, rowInfo, method);
      }
    }
  }

  nMoves++;
}


static void
MoveColToLeft(char **xy, int y, int startFunc, int lastFunc,
              int startVar, int lastVar, RcInfo_t *rowInfo,
              RcInfo_t *colInfo, int *rowOrder, int *colOrder,
              RcListInfo_t *candList, int method)
{
  int   i, j, s, t;
  int   row, col, tmpCol;

  assert(y >= startVar && y <= lastVar);

  col = colOrder[y];
  if (y != startVar) {
    /* Change col order */
    for (j = y; j > startVar; j--) {
      colOrder[j] = colOrder[j - 1];
      colInfo[colOrder[j]].pos = j;
    }
    colOrder[startVar] = col;
    colInfo[col].pos = startVar;
  }

  s = rowInfo[colInfo[col].gMin].pos;
  t = rowInfo[colInfo[col].gMax].pos;
  for (i = s; i <= t; i++) {
    row = rowOrder[i];
    if (!xy[row][col])
      continue;
    rowInfo[row].prevG = startVar;
    if (rowInfo[row].lNum == 1) {
      rowInfo[row].lMin = col;
      rowInfo[row].lMax = col;
    } else {
      if (col == rowInfo[row].lMin) {
        for (j = y + 1; j <= colInfo[rowInfo[row].lMax].pos; j++) {
          tmpCol = colOrder[j];
          if (xy[row][tmpCol]) {
            rowInfo[row].lMin = tmpCol;
            break;
          }
        }
      } else if (col == rowInfo[row].lMax) {
        for (j = y; j >= colInfo[rowInfo[row].lMin].pos; j--) {
          tmpCol = colOrder[j];
          if (xy[row][tmpCol]) {
            rowInfo[row].lMax = tmpCol;
            break;
          }
        }
      }
    }

    if (rowInfo[row].gNum == 1) {
      rowInfo[row].gMin = col;
      rowInfo[row].gMax = col;
    } else {
      if (col == rowInfo[row].gMin)
        rowInfo[row].gMin = col;
      else if (col == rowInfo[row].gMax) {
        rowInfo[row].gMax = rowInfo[row].lMax;
        if (startVar <= colInfo[rowInfo[row].gMin].pos)
          rowInfo[row].gMin = col;
      } else if (startVar <= colInfo[rowInfo[row].gMin].pos)
        rowInfo[row].gMin = col;
    }

    if (colInfo[col].lNum) {
      if (i == colInfo[col].prevG)
        i = rowInfo[colInfo[col].lMin].pos - 1;
      else if (i == rowInfo[colInfo[col].lMax].pos)
        i = colInfo[col].nextG - 1;
    }
  }

  s = rowInfo[colInfo[col].lMin].pos;
  t = rowInfo[colInfo[col].lMax].pos;
  for (i = s; i <= t; i++) {
    row = rowOrder[i];
    if (xy[row][col]) {
      rowInfo[row].lNum--;
      if (candList) {
        if (rowInfo[row].lNum == 0)
          SortedListDelete(candList, row);
        else if (rowInfo[row].lNum == 1)
          SortedListInsert(candList, row, col, colInfo, method);
      }
    }
  }

  nMoves++;
}


static void
MoveRowToBottom(char **xy, int x, int startFunc, int lastFunc,
                int startVar, int lastVar, RcInfo_t *rowInfo,
                RcInfo_t *colInfo, int *rowOrder, int *colOrder,
                RcListInfo_t *candList, int method)
{
  int   i, j, s, t;
  int   row, col, tmpRow;

  assert(x >= startFunc && x <= lastFunc);

  row = rowOrder[x];
  if (x != lastFunc) {
    /* Change row order */
    for (i = x; i < lastFunc; i++) {
      rowOrder[i] = rowOrder[i + 1];
      rowInfo[rowOrder[i]].pos = i;
    }
    rowOrder[lastFunc] = row;
    rowInfo[row].pos = lastFunc;
  }

  s = colInfo[rowInfo[row].gMin].pos;
  t = colInfo[rowInfo[row].gMax].pos;
  for (j = s; j <= t; j++) {
    col = colOrder[j];
    if (!xy[row][col])
      continue;
    colInfo[col].nextG = lastFunc;
    if (colInfo[col].lNum == 1) {
      colInfo[col].lMin = row;
      colInfo[col].lMax = row;
    } else {
      if (row == colInfo[col].lMin) {
        for (i = x; i <= rowInfo[colInfo[col].lMax].pos; i++) {
          tmpRow = rowOrder[i];
          if (xy[tmpRow][col]) {
            colInfo[col].lMin = tmpRow;
            break;
          }
        }
      } else if (row == colInfo[col].lMax) {
        for (i = x - 1; i >= rowInfo[colInfo[col].lMin].pos; i--) {
          tmpRow = rowOrder[i];
          if (xy[tmpRow][col]) {
            colInfo[col].lMax = tmpRow;
            break;
          }
        }
      }
    }

    if (colInfo[col].gNum == 1) {
      colInfo[col].gMin = row;
      colInfo[col].gMax = row;
    } else {
      if (row == colInfo[col].gMax)
        colInfo[col].gMax = row;
      else if (row == colInfo[col].gMin) {
        colInfo[col].gMin = colInfo[col].lMin;
        if (lastFunc >= rowInfo[colInfo[col].gMax].pos)
          colInfo[col].gMax = row;
      } else if (lastFunc >= rowInfo[colInfo[col].gMax].pos)
        colInfo[col].gMax = row;
    }

    if (rowInfo[row].lNum) {
      if (j == rowInfo[row].prevG)
        j = colInfo[rowInfo[row].lMin].pos - 1;
      else if (j == colInfo[rowInfo[row].lMax].pos)
        j = rowInfo[row].nextG - 1;
    }
  }

  s = colInfo[rowInfo[row].lMin].pos;
  t = colInfo[rowInfo[row].lMax].pos;
  for (j = s; j <= t; j++) {
    col = colOrder[j];
    if (xy[row][col]) {
      colInfo[col].lNum--;
      if (candList) {
        if (colInfo[col].lNum == 0)
          SortedListDelete(candList, col);
        else if (colInfo[col].lNum == 1)
          SortedListInsert(candList, col, row, rowInfo, method);
      }
    }
  }

  nMoves++;
}


static void
MoveColToRight(char **xy, int y, int startFunc, int lastFunc,
               int startVar, int lastVar, RcInfo_t *rowInfo,
               RcInfo_t *colInfo, int *rowOrder, int *colOrder,
               RcListInfo_t *candList, int method)
{
  int   i, j, s, t;
  int   row, col, tmpCol;

  assert(y >= startVar && y <= lastVar);

  col = colOrder[y];
  if (y != lastVar) {
    /* Change col order */
    for (j = y; j < lastVar; j++) {
      colOrder[j] = colOrder[j + 1];
      colInfo[colOrder[j]].pos = j;
    }
    colOrder[lastVar] = col;
    colInfo[col].pos = lastVar;
  }

  s = rowInfo[colInfo[col].gMin].pos;
  t = rowInfo[colInfo[col].gMax].pos;
  for (i = s; i <= t; i++) {
    row = rowOrder[i];
    if (!xy[row][col])
      continue;
    rowInfo[row].nextG = lastVar;
    if (rowInfo[row].lNum == 1) {
      rowInfo[row].lMin = col;
      rowInfo[row].lMax = col;
    } else {
      if (col == rowInfo[row].lMin) {
        for (j = y; j <= colInfo[rowInfo[row].lMax].pos; j++) {
          tmpCol = colOrder[j];
          if (xy[row][tmpCol]) {
            rowInfo[row].lMin = tmpCol;
            break;
          }
        }
      } else if (col == rowInfo[row].lMax) {
        for (j = y - 1; j >= colInfo[rowInfo[row].lMin].pos; j--) {
          tmpCol = colOrder[j];
          if (xy[row][tmpCol]) {
            rowInfo[row].lMax = tmpCol;
            break;
          }
        }
      }
    }

    if (rowInfo[row].gNum == 1) {
      rowInfo[row].gMin = col;
      rowInfo[row].gMax = col;
    } else {
      if (col == rowInfo[row].gMax)
        rowInfo[row].gMax = col;
      else if (col == rowInfo[row].gMin) {
        rowInfo[row].gMin = rowInfo[row].lMin;
        if (lastVar >= colInfo[rowInfo[row].gMax].pos)
          rowInfo[row].gMax = col;
      } else if (lastVar >= colInfo[rowInfo[row].gMax].pos)
        rowInfo[row].gMax = col;
    }

    if (colInfo[col].lNum) {
      if (i == colInfo[col].prevG)
        i = rowInfo[colInfo[col].lMin].pos - 1;
      else if (i == rowInfo[colInfo[col].lMax].pos)
        i = colInfo[col].nextG - 1;
    }
  }

  s = rowInfo[colInfo[col].lMin].pos;
  t = rowInfo[colInfo[col].lMax].pos;
  for (i = s; i <= t; i++) {
    row = rowOrder[i];
    if (xy[row][col]) {
      rowInfo[row].lNum--;
      if (candList) {
        if (rowInfo[row].lNum == 0)
          SortedListDelete(candList, row);
        else if (rowInfo[row].lNum == 1)
          SortedListInsert(candList, row, col, colInfo, method);
      }
    }
  }

  nMoves++;
}


static void
PrintMatrix(char **xy, int nRows, int nCols, int *rowOrder, int *colOrder,
            RcInfo_t *rowInfo, RcInfo_t *colInfo,
            int startFunc, int lastFunc, int startVar, int lastVar)
{
  int   i, j, x, y;
  int   row, col, index;
  char  line[2048];
  bdd_t *nsVar;

  line[nCols] = '\0';
  for (x = 0; x < nRows; x++) {
    if (startFunc != lastFunc && startVar != lastVar &&
        x == startFunc && x != 0) {
      for (y = 0; y < nCols; y++) {
        if (y == startVar || y == lastVar)
          line[y] = '%';
        else
          line[y] = '=';
      }
      fprintf(vis_stdout, "%s\n", line);
    }
    i = rowOrder[x];
    for (y = 0; y < nCols; y++) {
      j = colOrder[y];
      if (xy[i][j])
        line[y] = '1';
      else
        line[y] = '.';
    }
    fprintf(vis_stdout, "%s\n", line);
    if (startFunc != lastFunc && startVar != lastVar &&
        x == lastFunc && x < nRows - 1) {
      for (y = 0; y < nCols; y++) {
        if (y == startVar || y == lastVar)
          line[y] = '%';
        else
          line[y] = '=';
      }
      fprintf(vis_stdout, "%s\n", line);
    }
  }
  fprintf(vis_stdout, "row order :");
  for (x = 0; x < nRows; x++) {
    row = rowOrder[x];
    fprintf(vis_stdout, " %d(", row);
    if (rowInfo[row].data.row.nsVarBddArray) {
      for (i = 0; i < array_n(rowInfo[row].data.row.nsVarBddArray); i++) {
        nsVar = array_fetch(bdd_t *, rowInfo[row].data.row.nsVarBddArray, i);
        index = (int)bdd_top_var_id(nsVar);
        if (i == 0)
          fprintf(vis_stdout, "%d", index);
        else
          fprintf(vis_stdout, ",%d", index);
      }
    }
    fprintf(vis_stdout, ")");
  }
  fprintf(vis_stdout, "\n");
  fprintf(vis_stdout, "col order :");
  for (y = 0; y < nCols; y++) {
    col = colOrder[y];
    index = (int)bdd_top_var_id(colInfo[col].data.col.var);
    fprintf(vis_stdout, " %d(%d)", col, index);
  }
  fprintf(vis_stdout, "\n");
  fprintf(vis_stdout, "#active rows : %d-%d\n", startFunc, lastFunc);
  fprintf(vis_stdout, "#active cols : %d-%d\n", startVar, lastVar);
}


static void
PrintMatrixWithCluster(char **xy, ClusterList_t *headCluster, int nCols,
                        int *rowOrder, int *colOrder,
                        Img_DirectionType direction)
{
  int           i, j, x, y;
  char          line[2048];
  ClusterList_t *list, *tailCluster;

  line[nCols] = '\0';
  if (direction == Img_Forward_c) {
    tailCluster = headCluster;
    list = headCluster;
    while (list->next) {
      tailCluster = list->next;
      list = tailCluster;
    }
    list = tailCluster;
    while (list) {
      for (x = list->start; x <= list->end; x++) {
        i = rowOrder[x];
        for (y = 0; y < nCols; y++) {
          j = colOrder[y];
          if (xy[i][j])
            line[y] = '1';
          else
            line[y] = '.';
        }
        fprintf(vis_stdout, "%s\n", line);
      }
      list = list->prev;
      if (list) {
        for (y = 0; y < nCols; y++)
          line[y] = '-';
        fprintf(vis_stdout, "%s\n", line);
      }
    }
  } else {
    list = headCluster;
    while (list) {
      for (x = list->start; x <= list->end; x++) {
        i = rowOrder[x];
        for (y = 0; y < nCols; y++) {
          j = colOrder[y];
          if (xy[i][j])
            line[y] = '1';
          else
            line[y] = '.';
        }
        fprintf(vis_stdout, "%s\n", line);
      }
      list = list->next;
      if (list) {
        for (y = 0; y < nCols; y++)
          line[y] = '-';
        fprintf(vis_stdout, "%s\n", line);
      }
    }
  }
}


static void
PrintClusterMatrix(ClusterList_t *headCluster, int nCols,
                   int *colOrder, Img_DirectionType direction)
{
  int           j, y;
  char          line[2048];
  ClusterList_t *list, *tailCluster;

  line[nCols] = '\0';
  if (direction == Img_Forward_c) {
    tailCluster = headCluster;
    list = headCluster;
    while (list->next) {
      tailCluster = list->next;
      list = tailCluster;
    }
    list = tailCluster;
    while (list) {
      for (y = 0; y < nCols; y++) {
        j = colOrder[y];
        if (list->supports[j])
          line[y] = '1';
        else
          line[y] = '.';
      }
      fprintf(vis_stdout, "%s\n", line);
      list = list->prev;
      if (list) {
        for (y = 0; y < nCols; y++)
          line[y] = '-';
        fprintf(vis_stdout, "%s\n", line);
      }
    }
  } else {
    list = headCluster;
    while (list) {
      for (y = 0; y < nCols; y++) {
        j = colOrder[y];
        if (list->supports[j])
          line[y] = '1';
        else
          line[y] = '.';
      }
      fprintf(vis_stdout, "%s\n", line);
      list = list->next;
      if (list) {
        for (y = 0; y < nCols; y++)
          line[y] = '-';
        fprintf(vis_stdout, "%s\n", line);
      }
    }
  }
}


static void
CheckMatrix(char **xy, SccList_t *sccList, int nRows, int nCols,
            RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder,
            int startFunc, int lastFunc, int startVar, int lastVar, int local)
{
  int   i, j, x, y;
  int   lMin, lMax, gMin, gMax, lNum, gNum, prevG, nextG;
  int   error = 0;

  for (x = 0; x < nRows; x++) {
    i = rowOrder[x];
    if (x != rowInfo[i].pos) {
      fprintf(vis_stderr, "Error : wrong rowOrder (%d != %d) in Row %d\n",
             x, rowInfo[i].pos, x);
      error = 1;
    }

    lMin = gMin = nCols;
    lMax = gMax = -1;
    lNum = gNum = 0;
    prevG = -1;
    nextG = nCols;

    for (y = 0; y < nCols; y++) {
      j = colOrder[y];
      if (xy[i][j]) {
        gNum++;
        if (y < gMin)
          gMin = y;
        if (y > gMax)
          gMax = y;
        if (local) {
          if (x >= startFunc && x <= lastFunc &&
              y >= startVar && y <= lastVar) {
            if (y < lMin)
              lMin = y;
            if (y > lMax)
              lMax = y;
            lNum++;
          }
          if (!sccList ||
              (x >= sccList->startFunc && x <= sccList->lastFunc)) {
            if (y < startVar)
              prevG = y;
            if (y > lastVar && nextG == nCols)
              nextG = y;
          }
        }
      }
    }
    if (gNum != rowInfo[i].gNum) {
      fprintf(vis_stderr, "Error : wrong gNum (%d != %d) in Row %d\n",
             gNum, rowInfo[i].gNum, x);
      error = 1;
    }
    if (gNum > 0) {
      if (gMin != colInfo[rowInfo[i].gMin].pos) {
        fprintf(vis_stderr, "Error : wrong gMin (%d != %d) in Row %d\n",
               gMin, colInfo[rowInfo[i].gMin].pos, x);
        error = 1;
      }
      if (gMax != colInfo[rowInfo[i].gMax].pos) {
        fprintf(vis_stderr, "Error : wrong gMax (%d != %d) in Row %d\n",
               gMax, colInfo[rowInfo[i].gMax].pos, x);
        error = 1;
      }
    }
    if (local) {
      if (x >= startFunc && x <= lastFunc &&
          lNum != rowInfo[i].lNum) {
        fprintf(vis_stderr, "Error : wrong lNum (%d != %d) in Row %d\n",
               lNum, rowInfo[i].lNum, x);
        error = 1;
      }
      if (x >= startFunc && x <= lastFunc && lMin < nCols &&
          lMin != colInfo[rowInfo[i].lMin].pos) {
        fprintf(vis_stderr, "Error : wrong lMin (%d != %d) in Row %d\n",
               lMin, colInfo[rowInfo[i].lMin].pos, x);
        error = 1;
      }
      if (x >= startFunc && x <= lastFunc && lMax > -1 &&
          lMax != colInfo[rowInfo[i].lMax].pos) {
        fprintf(vis_stderr, "Error : wrong lMax (%d != %d) in Row %d\n",
               lMax, colInfo[rowInfo[i].lMax].pos, x);
        error = 1;
      }
      if (!sccList || (x >= sccList->startFunc && x <= sccList->lastFunc)) {
        if (prevG != rowInfo[i].prevG) {
          fprintf(vis_stderr, "Error : wrong prevG (%d != %d) in Row %d\n",
                  prevG, rowInfo[i].prevG, x);
          error = 1;
        }
        if (nextG != rowInfo[i].nextG) {
          fprintf(vis_stderr, "Error : wrong nextG (%d != %d) in Row %d\n",
                  nextG, rowInfo[i].nextG, x);
          error = 1;
        }
      }
    }
  }

  for (y = 0; y < nCols; y++) {
    j = colOrder[y];
    if (y != colInfo[j].pos) {
      fprintf(vis_stderr, "Error : wrong colOrder (%d != %d) in Col %d\n",
             y, colInfo[y].pos, y);
      error = 1;
    }

    lMin = gMin = nRows;
    lMax = gMax = -1;
    lNum = gNum = 0;
    prevG = -1;
    nextG = nRows;

    for (x = 0; x < nRows; x++) {
      i = rowOrder[x];
      if (xy[i][j]) {
        gNum++;
        if (x < gMin)
          gMin = x;
        if (x > gMax)
          gMax = x;
        if (local) {
          if (x >= startFunc && x <= lastFunc &&
              y >= startVar && y <= lastVar) {
            if (x < lMin)
              lMin = x;
            if (x > lMax)
              lMax = x;
            lNum++;
          }
          if (!sccList ||
              (y >= sccList->startVar && y <= sccList->lastVar)) {
            if (x < startFunc)
              prevG = x;
            if (x > lastFunc && nextG == nRows)
              nextG = x;
          }
        }
      }
    }
    if (gNum != colInfo[j].gNum) {
      fprintf(vis_stderr, "Error : wrong gNum (%d != %d) in Col %d\n",
             gNum, colInfo[j].gNum, y);
      error = 1;
    }
    if (gNum > 0) {
      if (gMin != rowInfo[colInfo[j].gMin].pos) {
        fprintf(vis_stderr, "Error : wrong gMin (%d != %d) in Col %d\n",
               gMin, rowInfo[colInfo[j].gMin].pos, y);
        error = 1;
      }
      if (gMax != rowInfo[colInfo[j].gMax].pos) {
        fprintf(vis_stderr, "Error : wrong gMax (%d != %d) in Col %d\n",
               gMax, rowInfo[colInfo[j].gMax].pos, y);
        error = 1;
      }
    }
    if (local) {
      if (y >= startVar && y <= lastVar &&
          lNum != colInfo[j].lNum) {
        fprintf(vis_stderr, "Error : wrong lNum (%d != %d) in Col %d\n",
               lNum, colInfo[j].lNum, y);
        error = 1;
      }
      if (y >= startVar && y <= lastVar && lMin < nRows &&
          lMin != rowInfo[colInfo[j].lMin].pos) {
        fprintf(vis_stderr, "Error : wrong lMin (%d != %d) in Col %d\n",
               lMin, rowInfo[colInfo[j].lMin].pos, y);
        error = 1;
      }
      if (y >= startVar && y <= lastVar && lMax > -1 &&
          lMax != rowInfo[colInfo[j].lMax].pos) {
        fprintf(vis_stderr, "Error : wrong lMax (%d != %d) in Col %d\n",
               lMax, rowInfo[colInfo[j].lMax].pos, y);
        error = 1;
      }
      if (!sccList || (y >= sccList->startVar && y <= sccList->lastVar)) {
        if (prevG != colInfo[j].prevG) {
          fprintf(vis_stderr, "Error : wrong prevG (%d != %d) in Col %d\n",
                  prevG, colInfo[j].prevG, y);
          error = 1;
        }
        if (nextG != colInfo[j].nextG) {
          fprintf(vis_stderr, "Error : wrong nextG (%d != %d) in Col %d\n",
                  nextG, colInfo[j].nextG, y);
          error = 1;
        }
      }
    }
  }

  assert(!error);
}


static void
CheckCluster(ClusterList_t *headCluster, int nCols,
                RcInfo_t *colInfo, int *colOrder)
{
  int   j, y, n;
  int   gMin, gMax, gNum;
  int   error = 0;
  ClusterList_t *list;

  n = 0;
  list = headCluster;
  while (list) {
    gMin = nCols;
    gMax = -1;
    gNum = 0;

    for (y = 0; y < nCols; y++) {
      j = colOrder[y];
      if (list->supports[j]) {
        gNum++;
        if (y < gMin)
          gMin = y;
        if (y > gMax)
          gMax = y;
      }
    }
    if (gNum != list->nSupports) {
      fprintf(vis_stderr, "Error : wrong nSupports (%d != %d) in Cluster %d\n",
             gNum, list->nSupports, n);
      error = 1;
    }
    if (gNum > 0) {
      if (gMin != colInfo[list->minCol].pos) {
        fprintf(vis_stderr, "Error : wrong gMin (%d != %d) in Cluster %d\n",
               gMin, colInfo[list->minCol].pos, n);
        error = 1;
      }
      if (gMax != colInfo[list->maxCol].pos) {
        fprintf(vis_stderr, "Error : wrong gMax (%d != %d) in Cluster %d\n",
               gMax, colInfo[list->maxCol].pos, n);
        error = 1;
      }
    }
    n++;
    list = list->next;
  }

  assert(!error);
}


static void
WriteMatrix(FILE *fout, char **xy, int nRows, int nCols, int *rowOrder,
            int *colOrder, RcInfo_t *rowInfo, RcInfo_t *colInfo)
{
  int   x, y;
  int   row, col;

  for (x = 0; x < nRows; x++) {
    row = rowOrder[x];
    for (y = 0; y < nCols; y++) {
      col = colOrder[y];
      if (xy[row][col])
        fprintf(fout, "%d %d %d\n", y, x, colInfo[col].data.col.type);
    }
  }
}


static void
PrintRow(char **xy, int nRows, int nCols, int *rowOrder, int *colOrder,
         int from, int to)
{
  int   x, y, row, col;

  for (x = from; x <= to; x++) {
    row = rowOrder[x];
    fprintf(vis_stdout, "[%d]%d :", x, row);
    for (y = 0; y < nCols; y++) {
      col = colOrder[y];
      if (xy[row][col])
        fprintf(vis_stdout, " %d(%d)", y, col);
    }
    fprintf(vis_stdout, "\n");
  }
}


static void
PrintCol(char **xy, int nRows, int nCols, int *rowOrder, int *colOrder,
         int from, int to)
{
  int   x, y, row, col;

  for (y = from; y <= to; y++) {
    col = colOrder[y];
    fprintf(vis_stdout, "[%d]%d :", y, col);
    for (x = 0; x < nRows; x++) {
      row = rowOrder[x];
      if (xy[row][col])
        fprintf(vis_stdout, " %d(%d)", x, row);
    }
    fprintf(vis_stdout, "\n");
  }
}


static RcListInfo_t *
SortedListAlloc(void)
{
  RcListInfo_t  *listInfo;

  listInfo = ALLOC(RcListInfo_t, 1);
  memset(listInfo, 0, sizeof(RcListInfo_t));
  listInfo->table = st_init_table(st_numcmp, st_numhash);
  return(listInfo);
}


static void
SortedListFree(RcListInfo_t *candList)
{
  RcList_t      *candidate, *next;

  candidate = candList->head;
  while (candidate) {
    next = candidate->next;
    FREE(candidate);
    candidate = next;
  }
  st_free_table(candList->table);
  FREE(candList);
}


static void
SortedListInsert(RcListInfo_t *candList, int index,
                 int otherIndex, RcInfo_t *otherInfo, int method)
{
  RcList_t      *candidate, *cur;
  int           lNum, gNum, diffH, diffT;

  candidate = ALLOC(RcList_t, 1);
  candidate->index = index;
  candidate->otherIndex = otherIndex;
  st_insert(candList->table, (char *)(long)index, (char *)candidate);

  if (candList->num == 0) {
    candidate->next = NIL(RcList_t);
    candidate->prev = NIL(RcList_t);
    candList->cur = candidate;
    candList->head = candidate;
    candList->tail = candidate;
    candList->num = 1;
    candList->maxNum = 1;
    return;
  }

  candList->num++;
  if (candList->num > candList->maxNum)
    candList->maxNum = candList->num;

  if (method == 1) {
    diffH = abs(index - candList->head->index);
    diffT = abs(index - candList->tail->index);

    if (diffH < diffT) {
      cur = candList->head;
      while (cur) {
        if (index < cur->index) {
          candidate->next = cur;
          candidate->prev = cur->prev;
          if (cur->prev)
            cur->prev->next = candidate;
          cur->prev = candidate;
          if (cur == candList->head)
            candList->head = candidate;
          break;
        }
        cur = cur->next;
      }
      if (!cur) {
        candidate->next = NIL(RcList_t);
        candidate->prev = candList->tail;
        candList->tail->next = candidate;
        candList->tail = candidate;
      }
    } else {
      cur = candList->tail;
      while (cur) {
        if (index > cur->index) {
          candidate->next = cur->next;
          if (cur->next)
            cur->next->prev = candidate;
          candidate->prev = cur;
          cur->next = candidate;
          if (cur == candList->tail)
            candList->tail = candidate;
          break;
        }
        cur = cur->prev;
      }
      if (!cur) {
        candidate->next = candList->head;
        candidate->prev = NIL(RcList_t);
        candList->head->prev = candidate;
        candList->head = candidate;
      }
    }

    if (candList->cur) {
      if (candList->cur->index == candList->curIndex)
        return;
      else if (candList->cur->index > candList->curIndex) {
        cur = candList->cur->prev;
        while (cur) {
          if (cur->index <= candList->curIndex)
            break;
          candList->cur = cur;
          cur = cur->prev;
        }
      } else {
        if (candidate->index > candList->curIndex)
          candList->cur = candidate;
        else
          candList->cur = candList->head;
      }
    } else
      candList->cur = candList->head;
  } else if (method == 2) {
    lNum = otherInfo[otherIndex].lNum;
    gNum = otherInfo[otherIndex].gNum;
    diffH = abs(lNum - otherInfo[candList->head->otherIndex].lNum);
    diffT = abs(lNum - otherInfo[candList->tail->otherIndex].lNum);

    if (diffH < diffT) {
      cur = candList->head;
      while (cur) {
        if (lNum < otherInfo[cur->otherIndex].lNum ||
            (lNum == otherInfo[cur->otherIndex].lNum &&
             (gNum < otherInfo[cur->otherIndex].gNum ||
              (gNum == otherInfo[cur->otherIndex].gNum &&
                index < cur->index)))) {
          candidate->next = cur;
          candidate->prev = cur->prev;
          if (cur->prev)
            cur->prev->next = candidate;
          cur->prev = candidate;
          if (cur == candList->head)
            candList->head = candidate;
          break;
        }
        cur = cur->next;
      }
      if (!cur) {
        candidate->next = NIL(RcList_t);
        candidate->prev = candList->tail;
        candList->tail->next = candidate;
        candList->tail = candidate;
      }
    } else {
      cur = candList->tail;
      while (cur) {
        if (lNum > otherInfo[cur->otherIndex].lNum ||
            (lNum == otherInfo[cur->otherIndex].lNum &&
             (gNum > otherInfo[cur->otherIndex].gNum ||
              (gNum == otherInfo[cur->otherIndex].gNum &&
                index > cur->index)))) {
          candidate->next = cur->next;
          if (cur->next)
            cur->next->prev = candidate;
          candidate->prev = cur;
          cur->next = candidate;
          if (cur == candList->tail)
            candList->tail = candidate;
          break;
        }
        cur = cur->prev;
      }
      if (!cur) {
        candidate->next = candList->head;
        candidate->prev = NIL(RcList_t);
        candList->head->prev = candidate;
        candList->head = candidate;
      }
    }
  } else if (method == 3) {
    lNum = otherInfo[otherIndex].lNum;
    gNum = otherInfo[otherIndex].gNum;
    diffH = abs(lNum - otherInfo[candList->head->otherIndex].lNum);
    diffT = abs(lNum - otherInfo[candList->tail->otherIndex].lNum);

    if (diffH < diffT) {
      cur = candList->head;
      while (cur) {
        if (lNum < otherInfo[cur->otherIndex].lNum ||
            (lNum == otherInfo[cur->otherIndex].lNum &&
             (gNum > otherInfo[cur->otherIndex].gNum ||
              (gNum == otherInfo[cur->otherIndex].gNum &&
                otherInfo[index].pos < otherInfo[cur->index].pos)))) {
          candidate->next = cur;
          candidate->prev = cur->prev;
          if (cur->prev)
            cur->prev->next = candidate;
          cur->prev = candidate;
          if (cur == candList->head)
            candList->head = candidate;
          break;
        }
        cur = cur->next;
      }
      if (!cur) {
        candidate->next = NIL(RcList_t);
        candidate->prev = candList->tail;
        candList->tail->next = candidate;
        candList->tail = candidate;
      }
    } else {
      cur = candList->tail;
      while (cur) {
        if (lNum > otherInfo[cur->otherIndex].lNum ||
            (lNum == otherInfo[cur->otherIndex].lNum &&
             (gNum < otherInfo[cur->otherIndex].gNum ||
              (gNum == otherInfo[cur->otherIndex].gNum &&
                otherInfo[index].pos > otherInfo[cur->index].pos)))) {
          candidate->next = cur->next;
          if (cur->next)
            cur->next->prev = candidate;
          candidate->prev = cur;
          cur->next = candidate;
          if (cur == candList->tail)
            candList->tail = candidate;
          break;
        }
        cur = cur->prev;
      }
      if (!cur) {
        candidate->next = candList->head;
        candidate->prev = NIL(RcList_t);
        candList->head->prev = candidate;
        candList->head = candidate;
      }
    }
  } else if (method == 4) {
    lNum = otherInfo[index].lNum;
    cur = candList->tail;
    while (cur) {
      if (cur->otherIndex > 1 ||
          (cur->otherIndex == 1 && otherInfo[cur->index].lNum > lNum) ||
          (cur->otherIndex == 1 && otherInfo[cur->index].lNum == lNum &&
           otherInfo[index].pos > otherInfo[cur->index].pos)) {
        candidate->next = cur->next;
        if (cur->next)
          cur->next->prev = candidate;
        candidate->prev = cur;
        cur->next = candidate;
        if (cur == candList->tail)
          candList->tail = candidate;
        break;
      }
      cur = cur->prev;
    }
    if (!cur) {
      candidate->next = candList->head;
      candidate->prev = NIL(RcList_t);
      candList->head->prev = candidate;
      candList->head = candidate;
    }
  } else {
    fprintf(vis_stderr, "** mlp error: invalid sort mode %d\n", method);
  }

  candList->cur = candList->head;
}


static void
SortedListDelete(RcListInfo_t *candList, int index)
{
  RcList_t      *candidate;
  long          lindex = (long) index;

  if (st_lookup(candList->table, (char *)lindex, &candidate)) {
    if (candidate == candList->head)
      candList->head = candidate->next;
    else
      candidate->prev->next = candidate->next;
    if (candidate == candList->tail)
      candList->tail = candidate->prev;
    else
      candidate->next->prev = candidate->prev;
    if (candidate == candList->cur) {
      if (candidate->next)
        candList->cur = candidate->next;
      else
        candList->cur = candList->head;
    }
    st_delete(candList->table, &lindex, NULL);
    candList->num--;
    FREE(candidate);
  }
}


static void
SortedListMove(RcListInfo_t *candList, RcInfo_t *info, int index, int method)
{
  RcList_t      *cur, *prev, *left, *right;

  if (st_lookup(candList->table, (char *)(long)index, &cur)) {
    if (method == 1) {
    } else if (method == 2) {
    } else if (method == 3) {
      while (cur != candList->head) {
        prev = cur->prev;
        if (info[prev->index].lNum < info[cur->index].lNum ||
            (info[prev->index].lNum == info[cur->index].lNum &&
             (info[prev->index].gNum > info[cur->index].gNum ||
              (info[prev->index].gNum == info[cur->index].gNum &&
                info[prev->index].pos < info[cur->index].pos)))) {
          break;
        }

        /* swap */
        if (prev == candList->head)
          candList->head = cur;
        if (cur == candList->tail)
          candList->tail = prev;
        left = prev->prev;
        right = cur->next;
        if (left)
          left->next = cur;
        cur->next = prev;
        prev->next = right;
        if (right)
          right->prev = prev;
        prev->prev = cur;
        cur->prev = left;
      }
    } else if (method == 4) {
      while (cur != candList->head) {
        prev = cur->prev;
        if (prev->otherIndex > cur->otherIndex ||
            (prev->otherIndex == cur->otherIndex &&
             (info[prev->index].lNum > info[cur->index].lNum ||
              (info[prev->index].lNum == info[cur->index].lNum &&
                info[prev->index].pos < info[cur->index].pos)))) {
          break;
        }

        /* swap */
        if (prev == candList->head)
          candList->head = cur;
        if (cur == candList->tail)
          candList->tail = prev;
        left = prev->prev;
        right = cur->next;
        if (left)
          left->next = cur;
        cur->next = prev;
        prev->next = right;
        if (right)
          right->prev = prev;
        prev->prev = cur;
        cur->prev = left;
      }
    }
  }
}


static void
CheckSortedList(RcListInfo_t *candList, RcInfo_t *info, int method)
{
  RcList_t      *cur, *next, *candidate;
  int           error = 0;
  int           num = 0;

  cur = candList->head;
  while (cur) {
    num++;
    next = cur->next;
    if (next) {
      if (next->prev != cur) {
        error++;
        if (next->prev) {
          fprintf(vis_stderr, "Error : prev of %d is not %d (!= %d) in list\n",
                 next->index, cur->index, next->prev->index);
        } else {
          fprintf(vis_stderr, "Error : prev of %d is not %d (!= nil) in list\n",
                 next->index, cur->index);
        }
      }
      if (method == 1) {
      } else if (method == 2) {
      } else if (method == 3) {
        if (info[cur->index].lNum > info[next->index].lNum) {
          error++;
          fprintf(vis_stderr,
                  "Error : cur(%d) lNum(%d) > next(%d) lNum(%d) in row list\n",
                  cur->index, info[cur->index].lNum,
                  next->index, info[next->index].lNum);
        } else if (info[cur->index].lNum == info[next->index].lNum &&
                   info[cur->index].gNum < info[next->index].gNum) {
          error++;
          fprintf(vis_stderr,
                  "Error : cur(%d) gNum(%d) < next(%d) gNum(%d) in row list\n",
                  cur->index, info[cur->index].gNum,
                  next->index, info[next->index].gNum);
        } else if (info[cur->index].lNum == info[next->index].lNum &&
                   info[cur->index].gNum == info[next->index].gNum &&
                   info[cur->index].pos > info[next->index].pos) {
          error++;
          fprintf(vis_stderr,
                  "Error : cur(%d) pos(%d) > next(%d) pos(%d) in row list\n",
                  cur->index, info[cur->index].pos,
                  next->index, info[next->index].pos);
        }
      } else if (method == 4) {
        if (cur->otherIndex < next->otherIndex) {
          error++;
          fprintf(vis_stderr, 
              "Error : cur(%d) length(%d) < next(%d) length(%d) in col list\n",
                  cur->index, cur->otherIndex, next->index, next->otherIndex);
        } else if (cur->otherIndex == next->otherIndex &&
                   info[cur->index].lNum < info[next->index].lNum) {
          error++;
          fprintf(vis_stderr,
                  "Error : cur(%d) lNum(%d) < next(%d) lNum(%d) in col list\n",
                  cur->index, info[cur->index].lNum,
                  next->index, info[next->index].lNum);
        } else if (cur->otherIndex == next->otherIndex &&
                   info[cur->index].lNum == info[next->index].lNum &&
                   info[cur->index].pos > info[next->index].pos) {
          error++;
          fprintf(vis_stderr,
                  "Error : cur(%d) pos(%d) > next(%d) pos(%d) in col list\n",
                  cur->index, info[cur->index].pos,
                  next->index, info[next->index].pos);
        }
      }
    } else {
      if (cur != candList->tail) {
        error++;
        if (candList->tail) {
          fprintf(vis_stderr, "Error : different tail in list (%d != %d)\n",
                 cur->index, candList->tail->index);
        } else {
          fprintf(vis_stderr, "Error : different tail in list (%d != nil)\n",
                cur->index);
        }
      }
    }
    if (!st_lookup(candList->table, (char *)(long)cur->index, &candidate)) {
      error++;
      fprintf(vis_stderr, "Error : index %d not in table\n", cur->index);
    }
    cur = next;
  }
  if (num != candList->num) {
    error++;
    fprintf(vis_stderr,
            "Error : different number of elements in list (%d != %d)\n",
            num, candList->num);
  }
  if (num != st_count(candList->table)) {
    error++;
    fprintf(vis_stderr,
            "Error : different number of elements in table (%d != %d)\n",
            num, st_count(candList->table));
  }
  assert(!error);
}


static void
MlpCluster(mdd_manager *mddManager, char **xy, int nRows, int nCols,
                int nActiveRows, int nActiveCols,
                int *nClusterRows, int *nClusterCols,
                int *rowOrder, int *colOrder,
                RcInfo_t *rowInfo, RcInfo_t *colInfo,
                array_t *clusterArray, array_t *arraySmoothVarBddArray,
                Img_DirectionType direction, int *cRowOrder,
                array_t *nsVarBddArray, int *sccBorder,
                int *varPos, ImgTrmOption_t *option)
{
  int                   i, j, row, col;
  int                   ncRows, ncCols;
  int                   index, qVarPos;
  int                   s1, t1, s2;
  array_t               *smoothVarBddArray;
  mdd_t                 *nsVar, *relation, *cluster, *tempCluster;
  ClusterList_t         *list, *headCluster, *tailCluster, *nextList, *prevList;
  ClusterSortedList_t   *clusterSortedList;
  long                  initialTime, finalTime;
  int                   moveFlag = 1;
  array_t               *nonAppearingVarBddArray;
  char                  *existFlag;
  array_t               *supportArray, *tmpArray;
  int                   nVars = bdd_num_vars(mddManager);
  mdd_t                 *one;

  if (option->mlpVerbosity)
    initialTime = util_cpu_time();
  else
    initialTime = 0;

  ncCols = nActiveCols;
  headCluster = NIL(ClusterList_t);
  tailCluster = NIL(ClusterList_t);
  clusterSortedList = NIL(ClusterSortedList_t);
  if (option->mlpCluster == 0) {
    if (direction == Img_Forward_c) {
      for (i = nActiveRows - 1; i >= 0;) {
        row = rowOrder[i];

        list = ALLOC(ClusterList_t, 1);
        memset(list, 0, sizeof(ClusterList_t));
        list->flag = 3;
        list->start = i;
        list->end = i;
        list->supports = ALLOC(char, nCols);
        memcpy(list->supports, xy[row], sizeof(char) * nCols);
        list->minCol = rowInfo[row].gMin;
        list->maxCol = rowInfo[row].gMax;
        list->nSupports = MlpCountSupport(list, colOrder, ncCols);
        list->product = bdd_dup(rowInfo[row].data.row.func);
        list->nQuantifyVars = 0;
        list->affinity = 0.0;
        if (rowInfo[row].data.row.nsVarBddArray)
          list->nsVarBddArray = array_dup(rowInfo[row].data.row.nsVarBddArray);
        else
          list->nsVarBddArray = NIL(array_t);

        if (headCluster)
          tailCluster->next = list;
        else
          headCluster = list;
        list->next = NIL(ClusterList_t);
        list->prev = tailCluster;
        tailCluster = list;

        i--;

        if (sccBorder && sccBorder[i + 1] == 1)
          continue;
        if (bdd_size(list->product) >= option->clusterSize)
          continue;

        while (i >= 0) {
          row = rowOrder[i];
          relation = rowInfo[row].data.row.func;
          if (bdd_size(list->product) >= option->clusterSize)
            break;
          cluster = bdd_and(list->product, relation, 1, 1);
          if (bdd_size(cluster) <= option->clusterSize) {
            mdd_free(list->product);
            list->product = cluster;
            list->start = i;

            if (list->nsVarBddArray) {
              array_append(list->nsVarBddArray,
                rowInfo[row].data.row.nsVarBddArray);
            }

            list->nSupports = 0;
            list->minCol = nActiveCols;
            list->maxCol = -1;
            s1 = nActiveCols;
            t1 = -1;
            supportArray = mdd_get_bdd_support_ids(mddManager, list->product);
            for (j = 0; j < array_n(supportArray); j++) {
              index = array_fetch(int, supportArray, j);
              if (varPos[index] == 0) {
                if ((int)bdd_top_var_id(colInfo[0].data.col.var) != index)
                  continue;
              }
              index = varPos[index];
              list->supports[index] = 1;
              list->nSupports++;
              s2 = colInfo[index].pos;
              if (s2 < s1) {
                s1 = s2;
                list->minCol = index;
              }
              if (s2 > t1) {
                t1 = s2;
                list->maxCol = index;
              }
            }
            array_free(supportArray);

            i--;
            if (option->mlpDebug)
              CheckCluster(headCluster, ncCols, colInfo, colOrder);
          } else {
            mdd_free(cluster);
            break;
          }
        }

        ncCols = RemoveLocalVarsInCluster(mddManager, xy,
                                            list, nActiveRows, ncCols,
                                            rowInfo, colInfo,
                                            rowOrder, colOrder,
                                            moveFlag, option);
      }
    } else {
      for (i = 0; i < nActiveRows;) {
        row = rowOrder[i];

        list = ALLOC(ClusterList_t, 1);
        memset(list, 0, sizeof(ClusterList_t));
        list->flag = 3;
        list->start = i;
        list->end = i;
        list->supports = ALLOC(char, nCols);
        memcpy(list->supports, xy[row], sizeof(char) * nCols);
        list->minCol = rowInfo[row].gMin;
        list->maxCol = rowInfo[row].gMax;
        list->nSupports = MlpCountSupport(list, colOrder, ncCols);
        list->product = bdd_dup(rowInfo[row].data.row.func);
        list->nQuantifyVars = 0;
        list->affinity = 0.0;
        if (rowInfo[row].data.row.nsVarBddArray)
          list->nsVarBddArray = array_dup(rowInfo[row].data.row.nsVarBddArray);
        else
          list->nsVarBddArray = NIL(array_t);

        if (headCluster)
          tailCluster->next = list;
        else
          headCluster = list;
        list->next = NIL(ClusterList_t);
        list->prev = tailCluster;
        tailCluster = list;

        i++;

        if (sccBorder && sccBorder[i - 1] == 2)
          continue;
        if (bdd_size(list->product) >= option->clusterSize)
          continue;

        while (i < nActiveRows) {
          row = rowOrder[i];
          relation = rowInfo[row].data.row.func;
          if (bdd_size(list->product) >= option->clusterSize)
            break;
          cluster = bdd_and(list->product, relation, 1, 1);
          if (bdd_size(cluster) <= option->clusterSize) {
            mdd_free(list->product);
            list->product = cluster;
            list->end = i;

            if (list->nsVarBddArray) {
              array_append(list->nsVarBddArray,
                rowInfo[row].data.row.nsVarBddArray);
            }

            list->nSupports = 0;
            list->minCol = nActiveCols;
            list->maxCol = -1;
            s1 = nActiveCols;
            t1 = -1;
            supportArray = mdd_get_bdd_support_ids(mddManager, list->product);
            for (j = 0; j < array_n(supportArray); j++) {
              index = array_fetch(int, supportArray, j);
              if (varPos[index] == 0) {
                if ((int)bdd_top_var_id(colInfo[0].data.col.var) != index)
                  continue;
              }
              index = varPos[index];
              list->supports[index] = 1;
              list->nSupports++;
              s2 = colInfo[index].pos;
              if (s2 < s1) {
                s1 = s2;
                list->minCol = index;
              }
              if (s2 > t1) {
                t1 = s2;
                list->maxCol = index;
              }
            }
            array_free(supportArray);

            i++;
            if (option->mlpDebug)
              CheckCluster(headCluster, ncCols, colInfo, colOrder);
          } else {
            mdd_free(cluster);
            break;
          }
        }

        ncCols = RemoveLocalVarsInCluster(mddManager, xy,
                                            list, nActiveRows, ncCols,
                                            rowInfo, colInfo,
                                            rowOrder, colOrder,
                                            moveFlag, option);
      }
    }
  } else { /* option->mlpCluster == 1 */
    if (direction == Img_Forward_c) {
      for (i = 0; i < nActiveRows; i++) {
        row = rowOrder[i];
        list = ALLOC(ClusterList_t, 1);
        memset(list, 0, sizeof(ClusterList_t));

        if (headCluster && headCluster->flag != 3)
          list->flag = 1;
        else
          list->flag = 2;
        list->start = i;
        list->end = i;
        list->supports = ALLOC(char, nCols);
        memcpy(list->supports, xy[row], sizeof(char) * nCols);
        list->minCol = rowInfo[row].gMin;
        list->maxCol = rowInfo[row].gMax;
        list->nSupports = MlpCountSupport(list, colOrder, ncCols);
        list->product = bdd_dup(rowInfo[row].data.row.func);
        list->nQuantifyVars = MlpNumQuantifyVars(list, rowInfo, colInfo,
                                                 colOrder, ncCols);
        if (rowInfo[row].data.row.nsVarBddArray)
          list->nsVarBddArray = array_dup(rowInfo[row].data.row.nsVarBddArray);
        else
          list->nsVarBddArray = NIL(array_t);

        if (bdd_size(list->product) >= option->clusterSize) {
          if (headCluster) {
            ncCols = RecursiveCluster(mddManager, headCluster,
                                        clusterSortedList,
                                        xy, rowInfo, colInfo,
                                        rowOrder, colOrder, nActiveRows, ncCols,
                                        direction, varPos, moveFlag, option);
            if (option->mlpDebug)
              CheckCluster(headCluster, ncCols, colInfo, colOrder);
          }
          list->flag = 3;
          ncCols = RemoveLocalVarsInCluster(mddManager, xy,
                                            list, nActiveRows, ncCols,
                                            rowInfo, colInfo,
                                            rowOrder, colOrder,
                                            moveFlag, option);
          if (option->mlpDebug)
            CheckCluster(headCluster, ncCols, colInfo, colOrder);
          list->affinity = 0.0;
          clusterSortedList = NIL(ClusterSortedList_t);
        } else if (sccBorder && i > 0 && sccBorder[i] == 1) {
          ncCols = RecursiveCluster(mddManager, headCluster,
                                        clusterSortedList,
                                        xy, rowInfo, colInfo,
                                        rowOrder, colOrder, nActiveRows, ncCols,
                                        direction, varPos, moveFlag, option);
          if (option->mlpDebug)
            CheckCluster(headCluster, ncCols, colInfo, colOrder);
          list->flag = 2;
          list->affinity = 0.0;
          clusterSortedList = NIL(ClusterSortedList_t);
          if (option->mlpClusterSortedList) {
            clusterSortedList = ClusterSortedListInsert(clusterSortedList, list,
                                                option->mlpClusterQuantifyVars);
          }
        } else {
          if (headCluster && headCluster->flag != 3) {
            list->affinity = MlpSupportAffinity(list, headCluster, colInfo,
                                                colOrder, ncCols,
                                                option->mlpCluster);
          } else
            list->affinity = 0.0;
          if (option->mlpClusterSortedList) {
            clusterSortedList = ClusterSortedListInsert(clusterSortedList, list,
                                                option->mlpClusterQuantifyVars);
          }
        }

        list->next = headCluster;
        list->prev = NIL(ClusterList_t);
        if (headCluster) {
          if (headCluster->flag == 1)
            headCluster->flag = 0;
          headCluster->prev = list;
        }
        headCluster = list;
      }
    } else {
      for (i = nActiveRows - 1; i >= 0; i--) {
        row = rowOrder[i];
        list = ALLOC(ClusterList_t, 1);
        memset(list, 0, sizeof(ClusterList_t));

        if (headCluster && headCluster->flag != 3)
          list->flag = 1;
        else
          list->flag = 2;
        list->start = i;
        list->end = i;
        list->supports = ALLOC(char, nCols);
        memcpy(list->supports, xy[row], sizeof(char) * nCols);
        list->minCol = rowInfo[row].gMin;
        list->maxCol = rowInfo[row].gMax;
        list->nSupports = MlpCountSupport(list, colOrder, ncCols);
        list->product = bdd_dup(rowInfo[row].data.row.func);
        list->nQuantifyVars = MlpNumQuantifyVars(list, rowInfo, colInfo,
                                                 colOrder, ncCols);
        if (rowInfo[row].data.row.nsVarBddArray)
          list->nsVarBddArray = array_dup(rowInfo[row].data.row.nsVarBddArray);
        else
          list->nsVarBddArray = NIL(array_t);

        if ( bdd_size(list->product) >= option->clusterSize) {
          if (headCluster) {
            ncCols = RecursiveCluster(mddManager, headCluster,
                                        clusterSortedList,
                                        xy, rowInfo, colInfo,
                                        rowOrder, colOrder, nActiveRows, ncCols,
                                        direction, varPos, moveFlag, option);
            if (option->mlpDebug)
              CheckCluster(headCluster, ncCols, colInfo, colOrder);
          }
          list->flag = 3;
          ncCols = RemoveLocalVarsInCluster(mddManager, xy,
                                            list, nActiveRows, ncCols,
                                            rowInfo, colInfo,
                                            rowOrder, colOrder,
                                            moveFlag, option);
          if (option->mlpDebug)
            CheckCluster(headCluster, ncCols, colInfo, colOrder);
          list->affinity = 0.0;
          clusterSortedList = NIL(ClusterSortedList_t);
        } else if (sccBorder && i < (nActiveRows - 1) && sccBorder[i] == 2) {
          ncCols = RecursiveCluster(mddManager, headCluster,
                                        clusterSortedList,
                                        xy, rowInfo, colInfo,
                                        rowOrder, colOrder, nActiveRows, ncCols,
                                        direction, varPos, moveFlag, option);
          if (option->mlpDebug)
            CheckCluster(headCluster, ncCols, colInfo, colOrder);
          list->flag = 2;
          list->affinity = 0.0;
          clusterSortedList = NIL(ClusterSortedList_t);
          if (option->mlpClusterSortedList) {
            clusterSortedList = ClusterSortedListInsert(clusterSortedList, list,
                                                option->mlpClusterQuantifyVars);
          }
        } else {
          if (headCluster && headCluster->flag != 3) {
            list->affinity = MlpSupportAffinity(list, headCluster, colInfo,
                                                colOrder, ncCols,
                                                option->mlpCluster);
          } else
            list->affinity = 0.0;
          if (option->mlpClusterSortedList) {
            clusterSortedList = ClusterSortedListInsert(clusterSortedList, list,
                                                option->mlpClusterQuantifyVars);
          }
        }

        list->next = headCluster;
        list->prev = NIL(ClusterList_t);
        if (headCluster) {
          if (headCluster->flag == 1)
            headCluster->flag = 0;
          headCluster->prev = list;
        }
        headCluster = list;

        if (option->mlpClusterSortedList && option->mlpDebug) {
          int   n1, n2;

          n1 = CountClusterList(headCluster);
          n2 = CountClusterSortedList(clusterSortedList);
          assert(n1 == n2);
        }
      }
    }

    if (headCluster->flag == 1) {
      ncCols = RecursiveCluster(mddManager, headCluster, clusterSortedList,
                                xy, rowInfo, colInfo,
                                rowOrder, colOrder, nActiveRows, ncCols,
                                direction, varPos, moveFlag, option);
      if (option->mlpDebug)
        CheckCluster(headCluster, ncCols, colInfo, colOrder);
    }
  }

  list = headCluster;
  one = bdd_one(mddManager);
  while (list) {
    if (bdd_equal(list->product, one)) {
      assert(list->nSupports == 0);

      prevList = list->prev;
      nextList = list->next;
      if ((!prevList) && (!nextList)) /* unity relation */
        break;

      mdd_free(list->product);
      if (list->nsVarBddArray)
        array_free(list->nsVarBddArray);
      FREE(list->supports);

      if (nextList) {
        if (list->start < nextList->start)
          nextList->start = list->start;
        if (list->end > nextList->end)
          nextList->end = list->end;

        if (list == headCluster)
          headCluster = nextList;

        FREE(list);
      } else {
        if (list->start < prevList->start)
          prevList->start = list->start;
        if (list->end > prevList->end)
          prevList->end = list->end;
        FREE(list);
      }
      if (prevList)
        prevList->next = nextList;
      if (nextList)
        nextList->prev = prevList;
      list = nextList;
      continue;
    }
    list = list->next;
  }
  bdd_free(one);

  if (option->mlpVerbosity >= 3 && option->mlpPrintMatrix) {
    PrintMatrixWithCluster(xy, headCluster, ncCols, rowOrder, colOrder,
                           direction);
    fprintf(vis_stdout, "******** cluster matrix ************\n");
    PrintClusterMatrix(headCluster, ncCols, colOrder, direction);
  }
  if (option->mlpDebug)
    CheckCluster(headCluster, ncCols, colInfo, colOrder);

  ncRows = 0;
  nonAppearingVarBddArray = NIL(array_t);
  if ((direction == Img_Forward_c && option->mlpReorder) ||
      (direction == Img_Backward_c && option->mlpPreReorder) ||
      option->mlpPostProcess >= 2) {
    list = headCluster;
    while (list) {
      row = rowOrder[ncRows];

      mdd_free(rowInfo[row].data.row.func);
      rowInfo[row].data.row.func = list->product;
      list->product = NIL(mdd_t);

      if (xy[row])
        FREE(xy[row]);
      xy[row] = list->supports;
      list->supports = NIL(char);

      rowInfo[row].gNum = list->nSupports;
      rowInfo[row].gMin = list->minCol;
      rowInfo[row].gMax = list->maxCol;
      rowInfo[row].lNum = list->nSupports;
      rowInfo[row].lMin = list->minCol;
      rowInfo[row].lMax = list->maxCol;
      rowInfo[row].prevG = -1;
      rowInfo[row].nextG = ncCols;

      if (list->nsVarBddArray) {
        if (rowInfo[row].data.row.nsVarBddArray)
          array_free(rowInfo[row].data.row.nsVarBddArray);
        rowInfo[row].data.row.nsVarBddArray = list->nsVarBddArray;
        list->nsVarBddArray = NIL(array_t);
      }

      nextList = list->next;
      FREE(list);
      ncRows++;
      list = nextList;
    }

    if (direction == Img_Forward_c) {
      for (i = 0; i < ncRows; i++) {
        cRowOrder[i] = rowOrder[ncRows - i - 1];
        rowInfo[rowOrder[ncRows - i - 1]].pos = i;
      }
    } else {
      for (i = 0; i < ncRows; i++) {
        cRowOrder[i] = rowOrder[i];
        rowInfo[rowOrder[i]].pos = i;
      }
    }

    for (j = 0; j < ncCols; j++) {
      col = colOrder[j];
      colInfo[col].gNum = 0;
      colInfo[col].gMin = ncRows;
      colInfo[col].gMax = -1;
      for (i = 0; i < ncRows; i++) {
        row = cRowOrder[i];
        if (xy[row][col]) {
          colInfo[col].gNum++;
          if (colInfo[col].gMin == ncRows)
            colInfo[col].gMin = row;
          colInfo[col].gMax = row;
        }
      }
      colInfo[col].lNum = colInfo[col].gNum;
      colInfo[col].lMin = colInfo[col].gMin;
      colInfo[col].lMax = colInfo[col].gMax;
      colInfo[col].prevG = -1;
      colInfo[col].nextG = ncRows;
    }
  } else {
    if (direction == Img_Forward_c) {
      if (arraySmoothVarBddArray) {
        nonAppearingVarBddArray = array_fetch(array_t *,
                                              arraySmoothVarBddArray, 0);
        if (nCols > ncCols) {
          existFlag = ALLOC(char, nVars);
          memset(existFlag, 0, sizeof(char) * nVars);

          for (i = 0; i < array_n(nonAppearingVarBddArray); i++) {
            nsVar = array_fetch(bdd_t *, nonAppearingVarBddArray, i);
            index = (int)bdd_top_var_id(nsVar);
            existFlag[index] = 1;
          }

          for (i = ncCols; i < nCols; i++) {
            col = colOrder[i];
            index = (int)bdd_top_var_id(colInfo[col].data.col.var);
            if (!existFlag[index]) {
              array_insert_last(bdd_t *, nonAppearingVarBddArray,
                                bdd_dup(colInfo[col].data.col.var));
            }
          }

          FREE(existFlag);
        }
        qVarPos = ncCols - 1;
        if (qVarPos >= 0) {
          col = colOrder[qVarPos];
          while (colInfo[col].gNum == 0) {
            array_insert_last(bdd_t *, nonAppearingVarBddArray,
                                bdd_dup(colInfo[col].data.col.var));
            if (qVarPos == 0)
              break;
            qVarPos--;
            col = colOrder[qVarPos];
          }
        }
      } else
        qVarPos = 0; /* to avoid warning */
    } else {
      if (arraySmoothVarBddArray) {
        nonAppearingVarBddArray = array_fetch(array_t *,
                                                arraySmoothVarBddArray, 0);
        if (nCols > ncCols) {
          for (i = ncCols; i < nCols; i++) {
            col = colOrder[i];
            if (colInfo[col].data.col.type == 2) {
              array_insert_last(bdd_t *, nonAppearingVarBddArray,
                                bdd_dup(colInfo[col].data.col.var));
            }
          }
        }
        qVarPos = ncCols - 1;
        if (qVarPos >= 0) {
          col = colOrder[qVarPos];
          while (colInfo[col].gNum == 0) {
            if (colInfo[col].data.col.type == 2) {
              array_insert_last(bdd_t *, nonAppearingVarBddArray,
                                bdd_dup(colInfo[col].data.col.var));
            }
            if (qVarPos == 0)
              break;
            qVarPos--;
            col = colOrder[qVarPos];
          }
        }

        existFlag = ALLOC(char, nVars);
        memset(existFlag, 0, sizeof(char) * nVars);

        for (i = 0; i < array_n(nonAppearingVarBddArray); i++) {
          nsVar = array_fetch(bdd_t *, nonAppearingVarBddArray, i);
          index = (int)bdd_top_var_id(nsVar);
          existFlag[index] = 1;
        }

        if (nRows > nActiveRows) {
          for (i = nActiveRows; i < nRows; i++) {
            row = rowOrder[i];
            tmpArray = rowInfo[row].data.row.nsVarBddArray;
            for (j = 0; j < array_n(tmpArray); j++) {
              nsVar = array_fetch(bdd_t *, tmpArray, j);
              index = (int)bdd_top_var_id(nsVar);
              existFlag[index] = 1;
            }
          }
        }

        list = headCluster;
        while (list) {
          supportArray = mdd_get_bdd_support_ids(mddManager, list->product);
          for (i = 0; i < array_n(supportArray); i++) {
            index = array_fetch(int, supportArray, i);
            existFlag[index] = 1;
          }
          array_free(supportArray);
          list = list->next;
        }

        list = headCluster;
        while (list) {
          for (i = 0; i < array_n(list->nsVarBddArray); i++) {
            nsVar = array_fetch(bdd_t *, list->nsVarBddArray, i);
            index = (int)bdd_top_var_id(nsVar);
            if (!existFlag[index]) {
              tmpArray = array_alloc(mdd_t *, 0);
              for (j = 0; j < i; j++) {
                nsVar = array_fetch(bdd_t *, list->nsVarBddArray, j);
                array_insert_last(mdd_t *, tmpArray, nsVar);
              }
              for (j = i + 1; j < array_n(list->nsVarBddArray); j++) {
                nsVar = array_fetch(bdd_t *, list->nsVarBddArray, j);
                index = (int)bdd_top_var_id(nsVar);
                if (existFlag[index])
                  array_insert_last(mdd_t *, tmpArray, nsVar);
              }
              array_free(list->nsVarBddArray);
              list->nsVarBddArray = tmpArray;
              break;
            }
          }
          list = list->next;
        }

        for (i = 0; i < array_n(nsVarBddArray); i++) {
          nsVar = array_fetch(bdd_t *, nsVarBddArray, i);
          index = (int)bdd_top_var_id(nsVar);
          if (!existFlag[index])
            array_insert_last(mdd_t *, nonAppearingVarBddArray, bdd_dup(nsVar));
        }

        FREE(existFlag);
      } else
        qVarPos = 0; /* to avoid warning */
    }

    if (direction == Img_Backward_c && nRows > nActiveRows) {
      cluster = bdd_one(mddManager);
      if (arraySmoothVarBddArray)
        smoothVarBddArray = array_alloc(array_t *, 0);
      else
        smoothVarBddArray = NIL(array_t);
      for (i = nActiveRows; i < nRows; i++) {
        row = rowOrder[i];
        relation = rowInfo[row].data.row.func;
        if (bdd_is_tautology(relation, 1))
          continue;
        if (smoothVarBddArray) {
          tmpArray = rowInfo[row].data.row.nsVarBddArray;
          for (j = 0; j < array_n(tmpArray); j++) {
            nsVar = array_fetch(bdd_t *, tmpArray, j);
            array_insert_last(mdd_t *, smoothVarBddArray, bdd_dup(nsVar));
          }
        }
        tempCluster = bdd_and(cluster, relation, 1, 1);
        bdd_free(cluster);
        cluster = tempCluster;
      }
      if (bdd_is_tautology(cluster, 1)) {
        mdd_free(cluster);
        if (smoothVarBddArray)
          mdd_array_free(smoothVarBddArray);
      } else {
        array_insert_last(bdd_t *, clusterArray, cluster);
        if (arraySmoothVarBddArray) {
          array_insert_last(array_t *, arraySmoothVarBddArray,
                            smoothVarBddArray);
        }
        ncRows++;
      }
    }

    list = headCluster;
    while (list) {
      array_insert_last(mdd_t *, clusterArray, list->product);
      list->product = NIL(mdd_t);

      if (arraySmoothVarBddArray) {
        if (direction == Img_Forward_c) {
          smoothVarBddArray = array_alloc(array_t *, 0);
          if (qVarPos >= 0) {
            col = colOrder[qVarPos];
            while (rowInfo[colInfo[col].gMin].pos >= list->start &&
                   rowInfo[colInfo[col].gMin].pos <= list->end) {
              array_insert_last(bdd_t *, smoothVarBddArray,
                                bdd_dup(colInfo[col].data.col.var));
              if (qVarPos == 0)
                break;
              qVarPos--;
              col = colOrder[qVarPos];
            }
          }
          array_insert_last(array_t *, arraySmoothVarBddArray,
                            smoothVarBddArray);
        } else {
          smoothVarBddArray = array_alloc(array_t *, 0);
          for (i = 0; i < array_n(list->nsVarBddArray); i++) {
            nsVar = array_fetch(bdd_t *, list->nsVarBddArray, i);
            array_insert_last(bdd_t *, smoothVarBddArray, bdd_dup(nsVar));
          }

          if (list->nSupports > 0) {
            s1 = colInfo[list->minCol].pos;
            t1 = colInfo[list->maxCol].pos;
            for (i = s1; i <= t1; i++) {
              col = colOrder[i];
              if (colInfo[col].data.col.type == 2) {
                if (rowInfo[colInfo[col].gMax].pos >= list->start &&
                    rowInfo[colInfo[col].gMax].pos <= list->end) {
                  array_insert_last(bdd_t *, smoothVarBddArray,
                                    bdd_dup(colInfo[col].data.col.var));
                }
              }
            }
          }
          array_insert_last(array_t *, arraySmoothVarBddArray,
                            smoothVarBddArray);
        }
      }

      nextList = list->next;
      if (list->nsVarBddArray)
        array_free(list->nsVarBddArray);
      FREE(list->supports);
      FREE(list);
      ncRows++;
      list = nextList;
    }

    if (direction == Img_Forward_c && nRows > nActiveRows) {
      cluster = bdd_one(mddManager);
      for (i = nActiveRows; i < nRows; i++) {
        row = rowOrder[i];
        relation = rowInfo[row].data.row.func;
        if (bdd_is_tautology(relation, 1))
          continue;
        tempCluster = bdd_and(cluster, relation, 1, 1);
        bdd_free(cluster);
        cluster = tempCluster;
      }
      if (bdd_is_tautology(cluster, 1))
        mdd_free(cluster);
      else {
        array_insert_last(bdd_t *, clusterArray, cluster);
        if (arraySmoothVarBddArray) {
          smoothVarBddArray = array_alloc(array_t *, 0);
          array_insert_last(array_t *, arraySmoothVarBddArray,
                            smoothVarBddArray);
        }
        ncRows++;
      }
    }
  }

  *nClusterRows = ncRows;
  *nClusterCols = ncCols;

  if (option->mlpVerbosity) {
    finalTime = util_cpu_time();
    fprintf(vis_stdout, "time for clustering = %10g\n",
           (double)(finalTime - initialTime) / 1000.0);
  }
}


static int
MlpCountSupport(ClusterList_t *list, int *colOrder, int nActiveCols)
{
  int   i, col;
  int   nSupports = 0;

  for (i = 0; i < nActiveCols; i++) {
    col = colOrder[i];
    if (list->supports[col])
      nSupports++;
  }
  return(nSupports);
}


static float
MlpSupportAffinity(ClusterList_t *curList, ClusterList_t *nextList,
                RcInfo_t *colInfo, int *colOrder,
                int nActiveCols, int clusterMethod)
{
  ClusterList_t *minList, *maxList;
  int           i, col, s, t;
  int           nOverlaps = 0;
  float         affinity;

  if (curList->nSupports <= nextList->nSupports) {
    minList = curList;
    maxList = nextList;
  } else {
    minList = nextList;
    maxList = curList;
  }

  s = colInfo[minList->minCol].pos;
  t = colInfo[minList->maxCol].pos;
  for (i = s; i <= t; i++) {
    col = colOrder[i];
    if (minList->supports[col] && maxList->supports[col])
      nOverlaps++;
  }

  if (clusterMethod == 1)
    affinity = (float)nOverlaps / (float)minList->nSupports;
  else {
    affinity = (float)nOverlaps /
                (float)(minList->nSupports + maxList->nSupports - nOverlaps);
  }
  return(affinity);
}


static int
RecursiveCluster(mdd_manager *mddManager, ClusterList_t *headCluster,
                 ClusterSortedList_t *clusterSortedList,
                 char **xy, RcInfo_t *rowInfo, RcInfo_t *colInfo,
                 int *rowOrder, int *colOrder,
                 int nActiveRows, int nClusterCols,
                 Img_DirectionType direction, int *varPos,
                 int moveFlag, ImgTrmOption_t *option)
{
  ClusterSortedList_t   *sortedList, *nextSortedList, *prevSortedList;
  ClusterSortedList_t   *firstSortedList, *secondSortedList, *saveSortedList;
  ClusterSortedList_t   *tailSortedList1, *tailSortedList2;
  ClusterList_t         *list, *best;
  int                   i, n1, n2;
  int                   s1, t1, s2;
  mdd_t                 *product;
  int                   merge, index;
  array_t               *supportArray;

  if (headCluster->flag == 3)
    return(nClusterCols);

  assert(headCluster);
  if (option->mlpClusterSortedList && option->mlpDebug) {
    n1 = CountClusterList(headCluster);
    n2 = CountClusterSortedList(clusterSortedList);
    assert(n1 == n2);
  }

  if (headCluster->flag == 2) {
    assert(!headCluster->next || headCluster->next->flag == 3);
    headCluster->flag = 3;
    nClusterCols = RemoveLocalVarsInCluster(mddManager, xy, headCluster,
                                            nActiveRows, nClusterCols,
                                            rowInfo, colInfo,
                                            rowOrder, colOrder,
                                            moveFlag, option);
    if (clusterSortedList) {
      assert(!clusterSortedList->next);
      FREE(clusterSortedList);
    }
    return(nClusterCols);
  }

  merge = 0;
  best = NIL(ClusterList_t);
  if (option->mlpClusterMerge) {
    /* NEED */
  }

  if (!best) {
    if (clusterSortedList)
      best = clusterSortedList->list;
    else {
      best = headCluster;
      list = headCluster->next;
      while (list) {
        if (list->flag == 2)
          break;
        if (option->mlpClusterQuantifyVars) {
          if (list->nQuantifyVars < best->nQuantifyVars ||
              (list->nQuantifyVars == best->nQuantifyVars &&
               list->affinity > best->affinity)) {
            best = list;
          }
        } else {
          if (list->affinity < option->mlpAffinityThreshold) {
            list = list->next;
            continue;
          }
          if (list->affinity > best->affinity ||
              (list->affinity == best->affinity &&
               list->nQuantifyVars < best->nQuantifyVars)) {
            best = list;
          }
        }
        list = list->next;
      }
    }
  }
  assert(best->next && best->next->flag != 3);

  list = best->next;
  product = mdd_and(best->product, list->product, 1, 1);

  if (merge == 0 && bdd_size(product) > option->clusterSize) {
    mdd_free(product);

    firstSortedList = NIL(ClusterSortedList_t);
    secondSortedList = NIL(ClusterSortedList_t);
    if (option->mlpClusterSortedList) {
      tailSortedList1 = NIL(ClusterSortedList_t);
      tailSortedList2 = NIL(ClusterSortedList_t);
      sortedList = clusterSortedList;
      while (sortedList) {
        nextSortedList = sortedList->next;
        if ((best->start > list->start &&
             sortedList->list->start >= best->start) ||
            (best->start < list->start &&
             sortedList->list->start <= best->start)) {
          if (tailSortedList1)
            tailSortedList1->next = sortedList;
          else
            firstSortedList = sortedList;
          tailSortedList1 = sortedList;
        } else {
          if (tailSortedList2)
            tailSortedList2->next = sortedList;
          else
            secondSortedList = sortedList;
          tailSortedList2 = sortedList;
        }
        sortedList->next = NIL(ClusterSortedList_t);
        sortedList = nextSortedList;
      }
    }

    if (best == headCluster) {
      best->flag = 3;
      if (option->mlpClusterSortedList) {
        sortedList = firstSortedList;
        firstSortedList = firstSortedList->next;
        FREE(sortedList);
        assert(!firstSortedList);
      }
      nClusterCols = RemoveLocalVarsInCluster(mddManager, xy, best,
                                                nActiveRows, nClusterCols,
                                                rowInfo, colInfo,
                                                rowOrder, colOrder,
                                                moveFlag, option);
    } else {
      best->flag = 2;
      if (option->mlpClusterSortedList) {
        sortedList = firstSortedList;
        prevSortedList = NIL(ClusterSortedList_t);
        saveSortedList = NIL(ClusterSortedList_t);
        while (sortedList) {
          nextSortedList = sortedList->next;
          if (sortedList->list == best) {
            if (nextSortedList) {
              if (prevSortedList)
                prevSortedList->next = sortedList->next;
              else
                firstSortedList = sortedList->next;
            }
            saveSortedList = sortedList;
            sortedList->next = NIL(ClusterSortedList_t);
            sortedList = nextSortedList;
            continue;
          }
          prevSortedList = sortedList;
          sortedList = nextSortedList;
        }
        prevSortedList->next = saveSortedList;
      }
      nClusterCols = RecursiveCluster(mddManager, headCluster, firstSortedList,
                                        xy, rowInfo, colInfo,
                                        rowOrder, colOrder,
                                        nActiveRows, nClusterCols,
                                        direction, varPos, moveFlag, option);
    }
    if (option->mlpDebug) {
      CheckCluster(headCluster, nClusterCols, colInfo, colOrder);
      CheckMatrix(xy, NIL(SccList_t), nActiveRows, nClusterCols,
                    rowInfo, colInfo, rowOrder, colOrder,
                    0, nActiveRows - 1, 0, nClusterCols - 1, 0);
    }

    if (list->flag == 0) {
      list->flag = 1;
      nClusterCols = RecursiveCluster(mddManager, list, secondSortedList,
                                        xy, rowInfo, colInfo,
                                        rowOrder, colOrder,
                                        nActiveRows, nClusterCols,
                                        direction, varPos, moveFlag, option);
    } else {
      list->flag = 3;
      sortedList = secondSortedList;
      secondSortedList = secondSortedList->next;
      FREE(sortedList);
      assert(!secondSortedList);
      nClusterCols = RemoveLocalVarsInCluster(mddManager, xy, list,
                                                nActiveRows, nClusterCols,
                                                rowInfo, colInfo,
                                                rowOrder, colOrder,
                                                moveFlag, option);
    }
    if (option->mlpDebug) {
      CheckCluster(headCluster, nClusterCols, colInfo, colOrder);
      CheckMatrix(xy, NIL(SccList_t), nActiveRows, nClusterCols,
                    rowInfo, colInfo, rowOrder, colOrder,
                    0, nActiveRows - 1, 0, nClusterCols - 1, 0);
    }
  } else {
    if (merge == 0 && option->mlpClusterSortedList) {
      sortedList = clusterSortedList->next;
      FREE(clusterSortedList);
      clusterSortedList = sortedList;
    }

    mdd_free(best->product);
    best->product = product;

    if (list->flag == 2) {
      if (best->flag == 1)
        best->flag = 3;
      else
        best->flag = 2;
    }

    best->nSupports = 0;
    best->minCol = nClusterCols;
    best->maxCol = -1;
    s1 = nClusterCols;
    t1 = -1;
    supportArray = mdd_get_bdd_support_ids(mddManager, best->product);
    for (i = 0; i < array_n(supportArray); i++) {
      index = array_fetch(int, supportArray, i);
      if (varPos[index] == 0) {
        if ((int)bdd_top_var_id(colInfo[0].data.col.var) != index)
          continue;
      }
      index = varPos[index];
      best->supports[index] = 1;
      best->nSupports++;
      s2 = colInfo[index].pos;
      if (s2 < s1) {
        s1 = s2;
        best->minCol = index;
      }
      if (s2 > t1) {
        t1 = s2;
        best->maxCol = index;
      }
    }
    array_free(supportArray);

    if (best->nsVarBddArray)
      array_append(best->nsVarBddArray, list->nsVarBddArray);

    if (list->start < best->start)
      best->start = list->start;
    if (list->end > best->end)
      best->end = list->end;

    if (option->mlpClusterSortedList) {
      sortedList = clusterSortedList;
      prevSortedList = NIL(ClusterSortedList_t);
      if (merge == 0) {
        while (sortedList) {
          nextSortedList = sortedList->next;
          if (option->mlpClusterDynamic) {
            if (sortedList->list == list || sortedList->list == best->prev) {
              if (prevSortedList)
                prevSortedList->next = sortedList->next;
              else
                clusterSortedList = sortedList->next;
              FREE(sortedList);
              sortedList = nextSortedList;
              continue;
            }
          } else {
            if (sortedList->list == list) {
              best->affinity = list->affinity;
              sortedList->list = best;
              break;
            }
          }
          prevSortedList = sortedList;
          sortedList = nextSortedList;
        }
      } else {
        while (sortedList) {
          nextSortedList = sortedList->next;
          if (option->mlpClusterDynamic) {
            if (sortedList->list == best ||
                sortedList->list == list || sortedList->list == best->prev) {
              if (prevSortedList)
                prevSortedList->next = sortedList->next;
              else
                clusterSortedList = sortedList->next;
              FREE(sortedList);
              sortedList = nextSortedList;
              continue;
            }
          } else {
            if (sortedList->list == best) {
              if (prevSortedList)
                prevSortedList->next = sortedList->next;
              else
                clusterSortedList = sortedList->next;
              FREE(sortedList);
              sortedList = nextSortedList;
              continue;
            } else if (sortedList->list == list) {
              best->affinity = list->affinity;
              sortedList->list = best;
              break;
            }
          }
          prevSortedList = sortedList;
          sortedList = nextSortedList;
        }
      }
    }

    best->next = list->next;
    if (list->next)
      list->next->prev = best;

    mdd_free(list->product);
    if (list->nsVarBddArray)
      array_free(list->nsVarBddArray);
    FREE(list->supports);
    FREE(list);

    if (merge && best->flag != 3) {
      nClusterCols = RemoveLocalVarsInCluster(mddManager, xy, best,
                                                nActiveRows, nClusterCols,
                                                rowInfo, colInfo,
                                                rowOrder, colOrder,
                                                moveFlag, option);
    }

    if (best->flag == 3) {
      nClusterCols = RemoveLocalVarsInCluster(mddManager, xy, best,
                                                nActiveRows, nClusterCols,
                                                rowInfo, colInfo,
                                                rowOrder, colOrder,
                                                moveFlag, option);
    } else {
      if (option->mlpClusterDynamic) {
        if (best->flag == 2)
          best->affinity = 0.0;
        else {
          best->affinity = MlpSupportAffinity(best, best->next, colInfo,
                                                colOrder, nClusterCols,
                                                option->mlpCluster);
        }
        if (best != headCluster) {
          best->prev->affinity = MlpSupportAffinity(best->prev, best, colInfo,
                                                    colOrder, nClusterCols,
                                                    option->mlpCluster);
        }
      }
      if (option->mlpClusterQuantifyVars == 2) {
        best->nQuantifyVars = MlpNumQuantifyVars(best, rowInfo, colInfo,
                                                 colOrder, nClusterCols);
      }
      if (option->mlpClusterSortedList && option->mlpClusterDynamic) {
        clusterSortedList = ClusterSortedListInsert(clusterSortedList, best,
                                                option->mlpClusterQuantifyVars);
        if (best != headCluster) {
          clusterSortedList = ClusterSortedListInsert(clusterSortedList,
                                                best->prev,
                                                option->mlpClusterQuantifyVars);
        }
      }
    }

    if (headCluster->flag == 1) {
      if (option->mlpDebug) {
        CheckCluster(headCluster, nClusterCols, colInfo, colOrder);
        CheckMatrix(xy, NIL(SccList_t), nActiveRows, nClusterCols,
                    rowInfo, colInfo, rowOrder, colOrder,
                    0, nActiveRows - 1, 0, nClusterCols - 1, 0);
      }
      nClusterCols = RecursiveCluster(mddManager, headCluster,
                                        clusterSortedList,
                                        xy, rowInfo, colInfo,
                                        rowOrder, colOrder,
                                        nActiveRows, nClusterCols,
                                        direction, varPos, moveFlag, option);
    } else {
      if (!option->mlpClusterDynamic) {
        assert(!clusterSortedList->next);
        FREE(clusterSortedList);
      }
    }
  }

  return(nClusterCols);
}


static int
RemoveLocalVarsInCluster(mdd_manager *mddManager, char **xy,
                         ClusterList_t *list,
                         int nActiveRows, int nClusterCols,
                         RcInfo_t *rowInfo, RcInfo_t *colInfo,
                         int *rowOrder, int *colOrder,
                         int moveFlag, ImgTrmOption_t *option)
{
  int           i, j, k, row, col, otherCol;
  int           s1, t1, s2, t2, s3, t3;
  mdd_t         *product;
  array_t       *localVarBddArray = NIL(array_t);

  s1 = colInfo[list->minCol].pos;
  t1 = colInfo[list->maxCol].pos;
  for (i = t1; i >= s1; i--) {
    col = colOrder[i];
    if (colInfo[col].data.col.type == 2 &&
        rowInfo[colInfo[col].gMin].pos >= list->start &&
        rowInfo[colInfo[col].gMax].pos <= list->end) {
      if (!localVarBddArray)
        localVarBddArray = array_alloc(bdd_t *, 0);
      array_insert_last(bdd_t *, localVarBddArray, colInfo[col].data.col.var);

      list->nSupports--;
      if (list->nSupports == 0) {
        list->minCol = -1;
        list->maxCol = -1;
      } else if (list->nSupports == 1) {
        if (list->minCol == col)
          list->minCol = list->maxCol;
        else
          list->maxCol = list->minCol;
      } else if (list->nSupports > 1) {
        if (list->minCol == col) {
          for (j = i + 1; j <= t1; j++) {
            otherCol = colOrder[j];
            if (list->supports[otherCol]) {
              list->minCol = otherCol;
              break;
            }
          }
        } else if (list->maxCol == col) {
          for (j = i - 1; j >= s1; j--) {
            otherCol = colOrder[j];
            if (list->supports[otherCol]) {
              list->maxCol = otherCol;
              break;
            }
          }
        }
      }

      s2 = list->start;
      t2 = list->end;
      for (j = s2; j <= t2; j++) {
        row = rowOrder[j];
        if (xy[row][col]) {
          rowInfo[row].gNum--;
          if (rowInfo[row].gNum > 0) {
            s3 = colInfo[rowInfo[row].gMin].pos;
            t3 = colInfo[rowInfo[row].gMax].pos;
            if (rowInfo[row].gMin == col) {
              for (k = i + 1; k <= t3; k++) {
                otherCol = colOrder[k];
                if (xy[row][otherCol]) {
                  rowInfo[row].gMin = otherCol;
                  break;
                }
              }
            } else if (rowInfo[row].gMax == col) {
              for (k = i - 1; k >= s3; k--) {
                otherCol = colOrder[k];
                if (xy[row][otherCol]) {
                  rowInfo[row].gMax = otherCol;
                  break;
                }
              }
            }
          }
        }
      }

      for (j = i; j < nClusterCols - 1; j++) {
        colOrder[j] = colOrder[j + 1];
        colInfo[colOrder[j]].pos = j;
      }
      colOrder[nClusterCols - 1] = col;
      colInfo[col].pos = nClusterCols - 1;

      nClusterCols--;

      if (option->mlpDebug) {
        CheckCluster(list, nClusterCols, colInfo, colOrder);
        CheckMatrix(xy, NIL(SccList_t), nActiveRows, nClusterCols,
                    rowInfo, colInfo, rowOrder, colOrder,
                    0, nActiveRows - 1, 0, nClusterCols - 1, 0);
      }
    }
  }

  if (localVarBddArray) {
    product = bdd_smooth(list->product, localVarBddArray);
    mdd_free(list->product);
    list->product = product;
    array_free(localVarBddArray);
    UpdateDisapearingPsVarsInCluster(mddManager, xy,
                        nActiveRows, nClusterCols,
                        rowOrder, colOrder, rowInfo, colInfo,
                        list, moveFlag, option);
  }

  return(nClusterCols);
}


static int
MlpNumQuantifyVars(ClusterList_t *list, RcInfo_t *rowInfo, RcInfo_t *colInfo,
                   int *colOrder, int nClusterCols)
{
  int   i, s, t, col;
  int   nQuantifyVars = 0;

  s = colInfo[list->minCol].pos;
  t = colInfo[list->maxCol].pos;
  for (i = t; i >= s; i--) {
    col = colOrder[i];
    if (rowInfo[colInfo[col].gMin].pos >= list->start &&
        rowInfo[colInfo[col].gMin].pos <= list->end) {
      nQuantifyVars++;
    } else
      break;
  }

  return(nQuantifyVars);
}


static ClusterSortedList_t *
ClusterSortedListInsert(ClusterSortedList_t *clusterSortedList,
                        ClusterList_t *list, int useQuantifyVars)
{
  ClusterSortedList_t   *sortedList, *cur, *prev;

  sortedList = ALLOC(ClusterSortedList_t, 1);
  sortedList->list = list;

  if (!clusterSortedList) {
    sortedList->next = NIL(ClusterSortedList_t);
    return(sortedList);
  }

  cur = clusterSortedList;
  prev = NIL(ClusterSortedList_t);
  while (cur) {
    if (list->flag == 2) {
      prev = cur;
      cur = cur->next;
      continue;
    }
    if (useQuantifyVars) {
      if (cur->list->flag == 2 ||
          list->nQuantifyVars < cur->list->nQuantifyVars ||
          (list->nQuantifyVars == cur->list->nQuantifyVars &&
           list->affinity >= cur->list->affinity)) {
        sortedList->next = cur;
        if (cur == clusterSortedList)
          clusterSortedList = sortedList;
        else
          prev->next = sortedList;
        break;
      }
    } else {
      if (cur->list->flag == 2 ||
          list->affinity > cur->list->affinity ||
          (list->affinity == cur->list->affinity &&
           list->nQuantifyVars <= cur->list->nQuantifyVars)) {
        sortedList->next = cur;
        if (cur == clusterSortedList)
          clusterSortedList = sortedList;
        else
          prev->next = sortedList;
        break;
      }
    }
    prev = cur;
    cur = cur->next;
  }
  if (!cur) {
    prev->next = sortedList;
    sortedList->next = NIL(ClusterSortedList_t);
  }
  return(clusterSortedList);
}


static int
CountClusterList(ClusterList_t *clusterList)
{
  ClusterList_t *list;
  int           n = 0;

  list = clusterList;
  while (list) {
    if (list->flag == 3)
      break;
    n++;
    if (list->flag == 2)
      break;
    list = list->next;
  }

  return(n);
}


static int
CountClusterSortedList(ClusterSortedList_t *clusterSortedList)
{
  ClusterSortedList_t   *sortedList;
  int                   n = 0;

  sortedList = clusterSortedList;
  while (sortedList) {
    n++;
    sortedList = sortedList->next;
  }

  return(n);
}


static array_t *
CreateInitialCluster(mdd_manager *mddManager, array_t *relationArray,
                     ImgFunctionData_t *functionData, array_t *nsVarBddArray,
                     ImgTrmOption_t *option)
{
  array_t               *clusteredRelationArray;
  array_t               *partitionArray;
  Part_Subsystem_t      *partitionSubsystem;
  st_table              *vertexTable;
  int                   i, j;
  int                   mddId;
  long                  bddId;
  Ntk_Node_t            *latch;
  mdd_t                 *cluster, *product, *relation, *var;
  char                  *latchName;
  st_generator          *stGen;
  st_table              *id2relation;
  array_t               *bddIdArray;
  int                   nVars;
  char                  *nsVarFlag;
  array_t               *domainVars;

  if (array_n(functionData->domainVars) == array_n(functionData->rangeVars))
    domainVars = functionData->domainVars;
  else {
    domainVars = array_alloc(int, 0);
    for (i = 0; i < array_n(functionData->rangeVars); i++) {
      mddId = array_fetch(int, functionData->domainVars, i);
      array_insert_last(int, domainVars, mddId);
    }
  }
  partitionArray = Fsm_ArdcDecomposeStateSpace(functionData->network,
                                                domainVars,
                                                functionData->roots,
                                                NIL(Fsm_ArdcOptions_t));
  if (domainVars != functionData->domainVars)
    array_free(domainVars);

  nVars = bdd_num_vars(mddManager);
  nsVarFlag = ALLOC(char, nVars);
  memset(nsVarFlag, 0, sizeof(char) * nVars);
  for (i = 0; i < array_n(nsVarBddArray); i++) {
    var = array_fetch(mdd_t *, nsVarBddArray, i);
    bddId = (long) bdd_top_var_id(var);
    nsVarFlag[bddId] = 1;
  }
  clusteredRelationArray = array_alloc(mdd_t *, 0);
  id2relation = st_init_table(st_numcmp, st_numhash);
  for (i = 0; i < array_n(relationArray); i++) {
    relation = array_fetch(mdd_t *, relationArray, i);
    bddIdArray = mdd_get_bdd_support_ids(mddManager, relation);
    for (j = 0; j < array_n(bddIdArray); j++) {
      bddId = (long) array_fetch(int, bddIdArray, j);
      if (nsVarFlag[bddId]) {
        st_insert(id2relation, (char *)bddId, (char *)relation);
        break;
      }
    }

    /* the relation of intermediate variable */
    if (j == array_n(bddIdArray)) {
      array_insert_last(mdd_t *, clusteredRelationArray, mdd_dup(relation));
    }

    array_free(bddIdArray);
  }
  FREE(nsVarFlag);

  arrayForEachItem(Part_Subsystem_t *, partitionArray, i, partitionSubsystem) {
    cluster = mdd_one(mddManager);
    vertexTable = Part_PartitionSubsystemReadVertexTable(partitionSubsystem);
    st_foreach_item(vertexTable, stGen, &latchName, NIL(char *)) {
      latch = Ntk_NetworkFindNodeByName(functionData->network, latchName); 
      mddId = Ntk_NodeReadMddId(Ntk_NodeReadShadow(latch));
      bddIdArray = mdd_id_to_bdd_id_array(mddManager, mddId);
      for (j = 0; j < array_n(bddIdArray); j++) {
        bddId = (long) array_fetch(int, bddIdArray, j);
        if (st_lookup(id2relation, (char *)bddId, &relation)) {
          st_delete(id2relation, &bddId, NULL);
          product = mdd_and(cluster, relation, 1, 1);
          mdd_free(cluster);
          cluster = product;
        }
      }
      array_free(bddIdArray);
    }
    Part_PartitionSubsystemFree(partitionSubsystem);
    array_insert_last(mdd_t *, clusteredRelationArray, cluster);
  }
  array_free(partitionArray);

  st_foreach_item(id2relation, stGen, &bddId, &relation) {
    array_insert_last(mdd_t *, clusteredRelationArray, mdd_dup(relation));
  }
  st_free_table(id2relation);

/*
  if (option->mlpDebug) {
    mdd_t       *tr1, *tr2, *tmp;

    tr1 = mdd_one(mddManager);
    for (i = 0; i < array_n(relationArray); i++) {
      relation = array_fetch(mdd_t *, relationArray, i);
      tmp = mdd_and(tr1, relation, 1, 1);
      mdd_free(tr1);
      tr1 = tmp;
    }

    tr2 = mdd_one(mddManager);
    for (i = 0; i < array_n(clusteredRelationArray); i++) {
      relation = array_fetch(mdd_t *, clusteredRelationArray, i);
      tmp = mdd_and(tr2, relation, 1, 1);
      mdd_free(tr2);
      tr2 = tmp;
    }

    assert(mdd_equal(tr1, tr2));
    mdd_free(tr1);
    mdd_free(tr2);
  }
*/

  return(clusteredRelationArray);
}


static void
SortCol(char **xy, int nRows, int nCols,
        RcInfo_t *rowInfo, RcInfo_t *colInfo, int *rowOrder, int *colOrder)
{
  int   x, y, i, j, row, col, lastVar;
  int   otherRow, otherCol;

  lastVar = nCols - 1;

  for (x = nRows - 1; x >= 0; x--) {
    row = rowOrder[x];
    for (y = lastVar; y >= 0; y--) {
      col = colOrder[y];
      if (colInfo[col].gMin == row) {
        if (y == lastVar) {
          lastVar--;
          continue;
        }
        for (j = y; j < lastVar; j++) {
          colOrder[j] = colOrder[j + 1];
          colInfo[colOrder[j]].pos = j;
        }
        colOrder[lastVar] = col;
        colInfo[col].pos = lastVar;
        for (i = x; i < nRows; i++) {
          otherRow = rowOrder[i];
          if (colInfo[rowInfo[otherRow].gMax].pos < lastVar)
            rowInfo[otherRow].gMax = col;
          if (rowInfo[otherRow].gMin == col) {
            for (j = 0; j < nCols; j++) {
              otherCol = colOrder[j];
              if (xy[otherRow][otherCol]) {
                rowInfo[otherRow].gMin = otherCol;
                break;
              }
            }
          }
        }
        lastVar--;
      }
    }
  }
}


static void
UpdateDisapearingPsVars(mdd_manager *mddManager, char **xy,
                        int nActiveRows, int nActiveCols,
                        int *rowOrder, int *colOrder,
                        RcInfo_t *rowInfo, RcInfo_t *colInfo,
                        int row, ImgTrmOption_t *option)
{
  int           i, id, x, y, s, t, col, tmpRow;
  mdd_t         *relation;
  array_t       *supportArray;
  char          *supportIndex;
  int           nVars;

  nVars = bdd_num_vars(mddManager);
  relation = rowInfo[row].data.row.func;
  supportIndex = ALLOC(char, nVars);
  memset(supportIndex, 0, sizeof(char) * nVars);
  supportArray = mdd_get_bdd_support_ids(mddManager, relation);
  for (i = 0; i < array_n(supportArray); i++) {
    id = array_fetch(int, supportArray, i);
    supportIndex[id] = 1;
  }
  array_free(supportArray);

  s = colInfo[rowInfo[row].gMin].pos;
  t = colInfo[rowInfo[row].gMax].pos;
  x = rowInfo[row].pos;
  rowInfo[row].gNum = 0;
  rowInfo[row].gMin = nActiveCols;
  rowInfo[row].gMax = -1;
  for (y = s; y <= t; y++) {
    col = colOrder[y];
    if (!xy[row][col])
      continue;
    id = (int)bdd_top_var_id(colInfo[col].data.col.var);
    if (supportIndex[id]) {
      rowInfo[row].gNum++;
      if (rowInfo[row].gMin == nActiveCols)
        rowInfo[row].gMin = col;
      rowInfo[row].gMax = col;
    } else {
      xy[row][col] = 0;
      if (colInfo[col].gNum == 1) {
        colInfo[col].gNum = 0;
        colInfo[col].gMin = nActiveRows;
        colInfo[col].gMax = -1;
      } else if (colInfo[col].gNum == 2) {
        colInfo[col].gNum = 1;
        if (colInfo[col].gMin == row)
          colInfo[col].gMin = colInfo[col].gMax;
        else
          colInfo[col].gMax = colInfo[col].gMin;
      } else {
        colInfo[col].gNum--;
        if (row == colInfo[col].gMin) {
          for (i = x + 1; i <= rowInfo[colInfo[col].gMax].pos; i++) {
            tmpRow = rowOrder[i];
            if (xy[tmpRow][col]) {
              colInfo[col].gMin = tmpRow;
              break;
            }
          }
        } else if (row == colInfo[col].gMax) {
          for (i = x - 1; i >= rowInfo[colInfo[col].gMin].pos; i--) {
            tmpRow = rowOrder[i];
            if (xy[tmpRow][col]) {
              colInfo[col].gMax = tmpRow;
              break;
            }
          }
        }
      }
    }
  }

  FREE(supportIndex);
}


static void
UpdateDisapearingPsVarsInCluster(mdd_manager *mddManager, char **xy,
                        int nActiveRows, int nActiveCols,
                        int *rowOrder, int *colOrder,
                        RcInfo_t *rowInfo, RcInfo_t *colInfo,
                        ClusterList_t *list, int moveFlag,
                        ImgTrmOption_t *option)
{
  int           i, j, k, y, id, row, col, otherRow, otherCol;
  int           s1, t1, s2, t2, s3, t3, t4;
  mdd_t         *relation;
  array_t       *supportArray;
  char          *supportIndex;
  int           nVars;
  ClusterList_t *otherList;

  if (list->nSupports == 0)
    return;

  nVars = bdd_num_vars(mddManager);
  relation = list->product;
  supportIndex = ALLOC(char, nVars);
  memset(supportIndex, 0, sizeof(char) * nVars);
  supportArray = mdd_get_bdd_support_ids(mddManager, relation);
  for (i = 0; i < array_n(supportArray); i++) {
    id = array_fetch(int, supportArray, i);
    supportIndex[id] = 1;
  }
  array_free(supportArray);

  s1 = colInfo[list->minCol].pos;
  t1 = colInfo[list->maxCol].pos;
  for (y = t1; y >= s1; y--) {
    col = colOrder[y];
    if (!list->supports[col])
      continue;
    id = (int)bdd_top_var_id(colInfo[col].data.col.var);
    if (!supportIndex[id]) {
      list->nSupports--;
      list->supports[col] = 0;
      if (list->nSupports == 0) {
        list->minCol = -1;
        list->maxCol = -1;
      } else if (list->nSupports == 1) {
        if (list->minCol == col)
          list->minCol = list->maxCol;
        else
          list->maxCol = list->minCol;
      } else if (list->nSupports > 1) {
        if (list->minCol == col) {
          for (j = y + 1; j <= t1; j++) {
            otherCol = colOrder[j];
            if (list->supports[otherCol]) {
              list->minCol = otherCol;
              break;
            }
          }
        } else if (list->maxCol == col) {
          for (j = y - 1; j >= s1; j--) {
            otherCol = colOrder[j];
            if (list->supports[otherCol]) {
              list->maxCol = otherCol;
              break;
            }
          }
        }
      }

      for (j = list->start; j <= list->end; j++) {
        row = rowOrder[j];
        if (xy[row][col]) {
          xy[row][col] = 0;
          colInfo[col].gNum--;
          rowInfo[row].gNum--;
          if (rowInfo[row].gNum > 0) {
            s2 = colInfo[rowInfo[row].gMin].pos;
            t2 = colInfo[rowInfo[row].gMax].pos;
            if (rowInfo[row].gMin == col) {
              for (k = y + 1; k <= t2; k++) {
                otherCol = colOrder[k];
                if (xy[row][otherCol]) {
                  rowInfo[row].gMin = otherCol;
                  break;
                }
              }
            } else if (rowInfo[row].gMax == col) {
              for (k = y - 1; k >= s2; k--) {
                otherCol = colOrder[k];
                if (xy[row][otherCol]) {
                  rowInfo[row].gMax = otherCol;
                  break;
                }
              }
            }
          } else {
            rowInfo[row].gMin = -1;
            rowInfo[row].gMax = nActiveCols;
          }
        }
      }

      if (colInfo[col].gNum == 0) {
        colInfo[col].lNum = 0;
        colInfo[col].gMin = nActiveRows;
        colInfo[col].lMin = nActiveRows;
        colInfo[col].gMax = -1;
        colInfo[col].lMax = -1;
      } else if (colInfo[col].gNum == 1) {
        colInfo[col].lNum = 1;
        if (rowInfo[colInfo[col].gMax].pos > list->end) {
          colInfo[col].gMin = colInfo[col].gMax;
          colInfo[col].lMin = colInfo[col].gMax;
        } else {
          colInfo[col].gMax = colInfo[col].gMin;
          colInfo[col].lMax = colInfo[col].gMin;
        }
      } else {
        colInfo[col].lNum = colInfo[col].gNum;
        if (rowInfo[colInfo[col].gMin].pos >= list->start &&
            rowInfo[colInfo[col].gMin].pos <= list->end) {
          for (i = list->end + 1; i <= rowInfo[colInfo[col].gMax].pos; i++) {
            otherRow = rowOrder[i];
            if (xy[otherRow][col]) {
              colInfo[col].gMin = otherRow;
              colInfo[col].lMin = otherRow;
              break;
            }
          }
        } else if (rowInfo[colInfo[col].gMax].pos >= list->start &&
                   rowInfo[colInfo[col].gMax].pos <= list->end) {
          for (i = list->start - 1; i >= rowInfo[colInfo[col].gMin].pos; i--) {
            otherRow = rowOrder[i];
            if (xy[otherRow][col]) {
              colInfo[col].gMax = otherRow;
              colInfo[col].lMax = otherRow;
              break;
            }
          }
        }
      }

      if (moveFlag) {
        if (colInfo[col].gNum == 0) {
          if (y < nActiveCols - 1) {
            for (j = y; j < nActiveCols - 1; j++) {
              colOrder[j] = colOrder[j + 1];
              colInfo[colOrder[j]].pos = j;
            }
            colOrder[nActiveCols - 1] = col;
            colInfo[col].pos = nActiveCols - 1;
          }
        } else if (rowInfo[colInfo[col].gMin].pos > list->end) {
          if (list->prev && list->prev->start > list->start) {
            /* Image */
            otherList = list->prev;
            while (!otherList->supports[col])
              otherList = otherList->prev;
          } else if (list->next && list->next->start > list->start) {
            /* Preimage */
            otherList = list->next;
            while (!otherList->supports[col])
              otherList = otherList->next;
          } else
            otherList = NIL(ClusterList_t);

          if (otherList) {
            t2 = colInfo[otherList->maxCol].pos;
            if (y != t2) {
              for (j = y; j < t2; j++) {
                colOrder[j] = colOrder[j + 1];
                colInfo[colOrder[j]].pos = j;
              }
              colOrder[t2] = col;
              colInfo[col].pos = t2;
              otherList->maxCol = col;
              if (otherList->minCol == col) {
                for (j = y; j < t2; j++) {
                  otherCol = colOrder[j];
                  if (otherList->supports[otherCol]) {
                    otherList->minCol = otherCol;
                    break;
                  }
                }
              }

              s3 = rowInfo[colInfo[col].gMin].pos;
              t3 = rowInfo[colInfo[col].gMax].pos;
              for (j = s3; j <= t3; j++) {
                otherRow = rowOrder[j];
                if (xy[otherRow][col]) {
                  if (rowInfo[otherRow].gNum > 1) {
                    t4 = colInfo[rowInfo[otherRow].gMax].pos;
                    if (t4 <= t2)
                      rowInfo[otherRow].gMax = col;
                    if (rowInfo[otherRow].gMin == col) {
                      for (k = y; k <= t4; k++) {
                        otherCol = colOrder[k];
                        if (xy[otherRow][otherCol]) {
                          rowInfo[otherRow].gMin = col;
                          break;
                        }
                      }
                    }
                  }
                }
              }
            }
          }
        }
      }

      if (option->mlpDebug) {
        CheckCluster(list, nActiveCols, colInfo, colOrder);
        CheckMatrix(xy, NIL(SccList_t), nActiveRows, nActiveCols,
                    rowInfo, colInfo, rowOrder, colOrder,
                    0, nActiveRows - 1, 0, nActiveCols - 1, 0);
      }
    }
  }

  FREE(supportIndex);
}


static void
UpdateNonappearingNsVars(mdd_manager *mddManager, array_t *nsVarBddArray,
                        int nRows, RcInfo_t *rowInfo, int *rowOrder,
                        array_t *nonAppearingVarBddArray)
{
  int           i, j, k, row, index, nVars;
  mdd_t         *nsVar, *relation;
  char          *existFlag;
  array_t       *supportArray, *tmpArray;

  nVars = bdd_num_vars(mddManager);
  existFlag = ALLOC(char, nVars);
  memset(existFlag, 0, sizeof(char) * nVars);

  for (i = 0; i < nRows; i++) {
    row = rowOrder[i];
    relation = rowInfo[row].data.row.func;
    supportArray = mdd_get_bdd_support_ids(mddManager, relation);
    for (j = 0; j < array_n(supportArray); j++) {
      index = array_fetch(int, supportArray, j);
      existFlag[index] = 1;
    }
    array_free(supportArray);
  }

  for (i = 0; i < nRows; i++) {
    row = rowOrder[i];
    for (j = 0; j < array_n(rowInfo[row].data.row.nsVarBddArray); j++) {
      nsVar = array_fetch(bdd_t *, rowInfo[row].data.row.nsVarBddArray, j);
      index = (int)bdd_top_var_id(nsVar);
      if (!existFlag[index]) {
        tmpArray = array_alloc(mdd_t *, 0);
        for (k = 0; k < j; k++) {
          nsVar = array_fetch(bdd_t *, rowInfo[row].data.row.nsVarBddArray, k);
          array_insert_last(mdd_t *, tmpArray, nsVar);
        }
        for (k = j + 1; k < array_n(rowInfo[row].data.row.nsVarBddArray); k++) {
          nsVar = array_fetch(bdd_t *, rowInfo[row].data.row.nsVarBddArray, k);
          index = (int)bdd_top_var_id(nsVar);
          if (existFlag[index])
            array_insert_last(mdd_t *, tmpArray, nsVar);
        }
        array_free(rowInfo[row].data.row.nsVarBddArray);
        rowInfo[row].data.row.nsVarBddArray = tmpArray;
        break;
      }
    }
  }

  for (i = 0; i < array_n(nsVarBddArray); i++) {
    nsVar = array_fetch(bdd_t *, nsVarBddArray, i);
    index = (int)bdd_top_var_id(nsVar);
    if (!existFlag[index])
      array_insert_last(mdd_t *, nonAppearingVarBddArray, bdd_dup(nsVar));
  }

  FREE(existFlag);
}


static void
WriteOrder(FILE *fout, int nCols, int *colOrder, RcInfo_t *colInfo)
{
  int   i, col;
  mdd_t *var;
  char  *name;

  for (i = 0; i < nCols; i++) {
    col = colOrder[i];
    var = colInfo[col].data.col.var;
    name = mdd_read_var_name(var);
    fprintf(fout, "%s\n", name);
  }
}