lp_MDO.c

/*
    Minimum matrix inverse fill-in modules - interface for lp_solve v5.0+
   ----------------------------------------------------------------------------------
    Author:        Kjell Eikland 
    Contact:       kjell.eikland@broadpark.no 
    License terms: LGPL.
    
    Requires:      string.h, colamd.h, lp_lib.h

    Release notes:
    v1.0    1 September 2003    Preprocessing routines for minimum fill-in column 
                                ordering for inverse factorization using the open 
                                source COLAMD library.  Suitable for the dense parts
                                of both the product form and LU factorization inverse 
                                methods.
    v1.1    1 July 2004         Renamed from lp_colamdMDO to lp_MDO.                                

   ---------------------------------------------------------------------------------- 
*/

#include <string.h>
#include "commonlib.h"
#include "lp_lib.h"
#include "colamd.h"
#include "lp_MDO.h"

#ifdef FORTIFY
# include "lp_fortify.h"
#endif

STATIC MYBOOL includeMDO(MYBOOL *usedpos, int item)
{
/*  Legend:   TRUE            => A basic slack variable already in the basis
              FALSE           => A column free for being pivoted in
              AUTOMATIC+TRUE  => A row-singleton user column pivoted into the basis
              AUTOMATIC+FALSE => A column-singleton user column pivoted into the basis */

  /* Handle case where we are processing all columns */
  if(usedpos == NULL)
    return( TRUE );
    
  else {
  /* Otherwise do the selective case */
    MYBOOL test = usedpos[item];
#if 1
    return( test != TRUE );
#else
    test = test & TRUE;
    return( test == FALSE );
#endif
  }
}

STATIC int prepareMDO(lprec *lp, MYBOOL *usedpos, int *colorder, int *data, int *rowmap)
/* This routine prepares data structures for colamd().  It is called twice, the first
   time to count applicable non-zero elements by column, and the second time to fill in 
   the row indexes of the non-zero values from the first call.  Note that the colamd() 
   row index base is 0 (which suits lp_solve fine). */
{
  int     i, ii, j, k, kk;
  int     nrows = lp->rows+1, ncols = colorder[0];
  int     offset = 0, Bnz = 0, Tnz;
  MYBOOL  dotally = (MYBOOL) (rowmap == NULL);
  MATrec  *mat = lp->matA;
  REAL    hold;
  REAL    *value;
  int     *rownr;

  if(dotally)
    data[0] = 0;

  Tnz = nrows - ncols;
  for(j = 1; j <= ncols; j++) {
    kk = colorder[j];

    /* Process slacks */
    if(kk <= lp->rows) {
      if(includeMDO(usedpos, kk)) {
        if(!dotally)
          data[Bnz] = rowmap[kk]+offset;
        Bnz++;
      }
      Tnz++;
    }
    /* Process user columns */
    else {
      k = kk - lp->rows;
      i = mat->col_end[k-1];
      ii= mat->col_end[k];
      Tnz += ii-i;
#ifdef Paranoia
      if(i >= ii)
        lp->report(lp, SEVERE, "prepareMDO: Encountered empty basic column %d\n", k);
#endif

      /* Detect if we need to do phase 1 adjustments of zero-valued OF variable */
      rownr = &COL_MAT_ROWNR(i);
      value = &COL_MAT_VALUE(i);
      hold = 0;
      if((*rownr > 0) && includeMDO(usedpos, 0) && modifyOF1(lp, kk, &hold, 1.0)) {
        if(!dotally)
          data[Bnz] = offset;
        Bnz++;
      }
      /* Loop over all NZ-variables */
      for(; i < ii; 
          i++, value += matValueStep, rownr += matRowColStep) {
        if(!includeMDO(usedpos, *rownr))
          continue;
        /* See if we need to change phase 1 OF value */
        if(*rownr == 0) {
          hold = *value;
          if(!modifyOF1(lp, kk, &hold, 1.0)) 
            continue;
        }
        /* Tally uneliminated constraint row values */
        if(!dotally)
          data[Bnz] = rowmap[*rownr]+offset;
        Bnz++;
      }
    }
    if(dotally)
      data[j] = Bnz;
  }
  return( Tnz );
}

STATIC MYBOOL verifyMDO(lprec *lp, int *col_end, int *row_nr, int rowmax, int colmax)
{
  int i, j, n, err = 0;

  for(i = 1; i <= colmax; i++) {
    n = 0;
    for(j = col_end[i-1]; (j < col_end[i]) && (err == 0); j++, n++) {
      if(row_nr[j] < 0 || row_nr[j] > rowmax)
        err = 1;
      if(n > 0 && row_nr[j] <= row_nr[j-1])
        err = 2;
      n++;
    }
  }
  if(err != 0)
    lp->report(lp, SEVERE, "verifyMDO: Invalid MDO input structure generated (error %d)\n", err);
  return( (MYBOOL) (err == 0) );
}

void *mdo_calloc(size_t size, size_t count)
{
  return ( calloc(size, count) );
}
void mdo_free(void *mem)
{
  free( mem );
}


int __WINAPI getMDO(lprec *lp, MYBOOL *usedpos, int *colorder, int *size, MYBOOL symmetric)
{
  int    error = FALSE;
  int    nrows = lp->rows+1, ncols = colorder[0];
  int    i, j, kk, n;
  int    *col_end, *row_map = NULL;
  int    Bnz, Blen, *Brows = NULL;
  int    stats[COLAMD_STATS];
  double knobs[COLAMD_KNOBS];

 /* Tally the non-zero counts of the unused columns/rows of the 
    basis matrix and store corresponding "net" starting positions */
  allocINT(lp, &col_end, ncols+1, FALSE);
  n = prepareMDO(lp, usedpos, colorder, col_end, NULL);
  Bnz = col_end[ncols];

 /* Check that we have unused basic columns, otherwise skip analysis */  
  if(ncols == 0 || Bnz == 0) 
    goto Transfer;

 /* Get net number of rows and fill mapper */
  allocINT(lp, &row_map, nrows, FALSE);
  nrows = 0;
  for(i = 0; i <= lp->rows; i++) {
    row_map[i] = i-nrows;
   /* Increment eliminated row counter if necessary */
    if(!includeMDO(usedpos, i)) 
      nrows++;
  }
  nrows = lp->rows+1 - nrows;

 /* Store row indeces of non-zero values in the basic columns */
  Blen = colamd_recommended(Bnz, nrows, ncols);
  allocINT(lp, &Brows, Blen, FALSE);
  prepareMDO(lp, usedpos, colorder, Brows, row_map);
#ifdef Paranoia
  verifyMDO(lp, col_end, Brows, nrows, ncols);
#endif

 /* Compute the MDO */
#if 1
  colamd_set_defaults(knobs);
  knobs [COLAMD_DENSE_ROW] = 0.2+0.2 ;    /* default changed for UMFPACK */
  knobs [COLAMD_DENSE_COL] = knobs [COLAMD_DENSE_ROW];    
  if(symmetric && (nrows == ncols)) {
    MEMCOPY(colorder, Brows, ncols + 1);
    error = !symamd(nrows, colorder, col_end, Brows, knobs, stats, mdo_calloc, mdo_free);
  }
  else
    error = !colamd(nrows, ncols, Blen, Brows, col_end, knobs, stats);
#else
  if(symmetric && (nrows == ncols)) {
    MEMCOPY(colorder, Brows, ncols + 1);
    error = !symamd(nrows, colorder, col_end, Brows, knobs, stats, mdo_calloc, mdo_free);
  }
  else
    error = !colamd(nrows, ncols, Blen, Brows, col_end, (double *) NULL, stats);
#endif

 /* Transfer the estimated optimal ordering, adjusting for index offsets */
Transfer:
  if(error) 
    error = stats[COLAMD_STATUS];
  else {
    MEMCOPY(Brows, colorder, ncols + 1);
    for(j = 0; j < ncols; j++) {
      kk = col_end[j];
      n = Brows[kk+1];
      colorder[j+1] = n;
    }
  }

  /* Free temporary vectors */
  FREE(col_end);
  if(row_map != NULL)
    FREE(row_map);
  if(Brows != NULL)
    FREE(Brows);

  if(size != NULL)
    *size = ncols;
  return( error );
}