generic/solver/slv_stdcalls.c

/*  ASCEND modelling environment
    Copyright (C) 1998, 2006 Carnegie Mellon University
    Copyright (C) 1996 Benjamin Andrew Allan

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2, or (at your option)
    any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place - Suite 330,
    Boston, MA 02111-1307, USA.
*//*
    by Benjamin Andrew Allan
    5/19/96
    Last in CVS: $Revision: 1.28 $ $Date: 1998/06/16 16:53:04 $ $Author: mthomas $
*/

#include <utilities/ascConfig.h>
#include <compiler/compiler.h>
#include <utilities/ascMalloc.h>
#include <general/list.h>
#include "mtx.h"
#include "slv_types.h"
#include "var.h"
#include "rel.h"
#include "discrete.h"
#include "conditional.h"
#include "logrel.h"
#include "relman.h"
#include "logrelman.h"
#include "bnd.h"
#include "slv_common.h"
#include "linsol.h"
#include "linsolqr.h"
#include "slv_client.h"

/* headers of registered clients */
#include "slv0.h"
#include "slv1.h"
#include "slv2.h"
#include "slv3.h"
#include "slv6.h"
#include "slv7.h"
#include "slv8.h"
#include "slv9.h"
#include "slv9a.h"

#include "model_reorder.h"
#include "slv_stdcalls.h"
#include <general/mathmacros.h>

#define KILL 0 /* deleteme code. compile old code if kill = 1 */
#define NEEDSTOBEDONE 0 /* indicates code/comments that are not yet ready */
#define USECODE 0  /* Code in good shape but not used currently */

#define MIMDEBUG 0 /* slv_std_make_incidence_mtx debugging */
#define RIDEBUG 0 /* reindex debugging */
#define SBPDEBUG 0 /* slv_block_partition_real debugging */
#define SLBPDEBUG 0 /* slv_log_block_partition debugging */

/* global to get around the mr header (for tear_subreorder) */
static
enum mtx_reorder_method g_blockmethod = mtx_UNKNOWN;

/*
    Here we play some hefty Jacobian reordering games.

    What we want to happen eventually is as follows:
    
      - Get the free & incident pattern for include relations.
      - Output-assign the Jacobian.
      - BLT permute the Jacobian. If underspecified, fake rows
        to make things appear square.
      - For all leading square blocks apply reordering (kirk, etc),
      - If trailing rectangular block, "apply clever kirkbased scheme not known."???
    
    At present, we aren't quite so clever. We:
    
      - Get the free & incident pattern for include relations.
      - Output-assign the Jacobian.
      - BLT permute the square output assigned region of the Jacobian.
      - For all leading square blocks apply reordering (kirk, etc)
      - Collect the block list as part of the master data structure.
      - Set sindices for rels, vars as current rows/cols in matrix so
        that the jacobian is 'naturally' preordered for all solvers.
    
    Solvers are still free to reorder their own matrices any way they like.
    It's probably a dumb idea, though.
*/
int slv_std_make_incidence_mtx(slv_system_t sys, mtx_matrix_t mtx,
                               var_filter_t *vf,rel_filter_t *rf)
{
#if MIMDEBUG
  FILE *fp;
#endif
  int32 r, rlen,ord;
  struct rel_relation **rp;

  if (sys==NULL || mtx == NULL || vf == NULL || rf == NULL) {
    ERROR_REPORTER_HERE(ASC_PROG_ERR,"called with null");
    return 1;
  }
  rp = slv_get_solvers_rel_list(sys);
  assert(rp!=NULL);
  rlen = slv_get_num_solvers_rels(sys);
  ord = MAX(rlen,slv_get_num_solvers_vars(sys));
  if (ord > mtx_order(mtx)) {
    ERROR_REPORTER_HERE(ASC_PROG_ERR,"undersized matrix");
    return 2;
  }
  for (r=0; r < rlen; r++) {
    if (rel_apply_filter(rp[r],rf)) {
      relman_get_incidence(rp[r],vf,mtx);
    }
  }
#if MIMDEBUG
  fp = fopen("/tmp/mim.plot","w+");
  if (fp !=NULL) {
    mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
    fclose(fp);
  }
#endif
  return 0;
}


/**
    Orders the solvers_var list of the system to match the permutation 
    on the given mtx. Does not change the data in mtx. 

    @return 0 on success, 1 on out-of-memory

    It is assumed that the org cols of mtx == var list position. Only the
    vars in range lo to hi of the var list are affected and
    only these vars should appear in the same org column range of
    the mtx. This should not be called on blocks less than 3x3.
*/
static int reindex_vars_from_mtx(slv_system_t sys, int32 lo, int32 hi,
                                 const mtx_matrix_t mtx)
{
  struct var_variable **vtmp, **vp;
  int32 c,v,vlen;

  if (lo >= hi +1) return 0; /* job too small */
  vp = slv_get_solvers_var_list(sys);
  vlen = slv_get_num_solvers_vars(sys);
  /* on vtmp we DONT have the terminating null */
  vtmp = ASC_NEW_ARRAY(struct var_variable *,vlen);
  if (vtmp == NULL) {
    return 1;
  }
  /* copy pointers to vtmp in order desired and copy back rather than sort */
  for (c=lo;c<=hi;c++) {
    v = mtx_col_to_org(mtx,c);
    vtmp[c] = vp[v];
  }
  /* copying back and re-sindexing */
  for (c=lo;c<=hi;c++) {
    vp[c] = vtmp[c];
    var_set_sindex(vp[c],c);
  }
  ascfree(vtmp);
  return 0;
}
/**
    Orders the solvers_rel list of the system to match the permutation 
    on the given mtx. Does not change the data in mtx. 

    @return 0 on success, 1 on out-of-memory

    It is assumed that the org rows of mtx == rel list position. Only the
    rels in range lo to hi of the rel list are affected and
    only these rels should appear in the same org row range of
    the input mtx. This should not be called on blocks less than 3x3.
*/
static int reindex_rels_from_mtx(slv_system_t sys, int32 lo, int32 hi,
                                 const mtx_matrix_t mtx)
{
  struct rel_relation **rtmp, **rp;
  int32 c,v,rlen;

  rp = slv_get_solvers_rel_list(sys);
  rlen = slv_get_num_solvers_rels(sys);
  /* on rtmp we DONT have the terminating null */
  rtmp = ASC_NEW_ARRAY(struct rel_relation*,rlen);
  if (rtmp == NULL) {
    return 1;
  }
  /* copy pointers to rtmp in order desired and copy back rather than sort.
   * do this only in the row range of interest. */
  for (c=lo;c<=hi;c++) {
    v = mtx_row_to_org(mtx,c);
#if RIDEBUG
    if(c!=v){
      CONSOLE_DEBUG("Old rel sindex (org) %d becoming sindex (cur) %d\n",v,c);
    }
#endif
    rtmp[c] = rp[v];
  }
  /* copying back and re-sindexing */
  for (c=lo;c<=hi;c++) {
    rp[c] = rtmp[c];
    rel_set_sindex(rp[c],c);
  }
  ascfree(rtmp);
  return 0;
}

/**
    Perform var and rel reordering to achieve block form.

    @param uppertrianguler if non-zero, BUT form. if 0, make BLT form.

    @callergraph
*/
int slv_block_partition_real(slv_system_t sys,int uppertriangular){
#if SBPDEBUG
  FILE *fp;
#endif
  struct rel_relation **rp;
  struct var_variable **vp;
  mtx_region_t  *newblocks;
  dof_t *d;
  int32 nrow,ncol;
  int32 ncolpfix, nrowpun, vlenmnv;
  mtx_matrix_t mtx;
  int32 order, rank;
  int32 c,len,vlen,r,rlen;
  var_filter_t vf;
  rel_filter_t rf;

  /* CONSOLE_DEBUG("..."); */

  rp = slv_get_solvers_rel_list(sys);
  vp = slv_get_solvers_var_list(sys);
  rlen = slv_get_num_solvers_rels(sys);
  vlen = slv_get_num_solvers_vars(sys);
  if (rlen ==0 || vlen == 0) return 1;
  order = MAX(rlen,vlen);

  rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
  rf.matchvalue = (REL_ACTIVE);
  vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
  vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);

  /* nrow plus unincluded rels and ineqs */
  nrowpun = slv_count_solvers_rels(sys,&rf);
  ncolpfix = slv_count_solvers_vars(sys,&vf); /* ncol plus fixed vars */

  vf.matchbits = (VAR_SVAR);
  vf.matchvalue = 0;
  vlenmnv = vlen - slv_count_solvers_vars(sys,&vf);  /* vlen minus non vars */

  rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
  rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
  vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED | VAR_ACTIVE);
  vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);

  nrow = slv_count_solvers_rels(sys,&rf);
  ncol = slv_count_solvers_vars(sys,&vf);

  mtx = mtx_create();
  mtx_set_order(mtx,order);

  if (slv_std_make_incidence_mtx(sys,mtx,&vf,&rf)) {
    ERROR_REPORTER_HERE(ASC_PROG_ERR,"failure in creating incidence matrix.");
    mtx_destroy(mtx);
    return 1;
  }

  /* CONSOLE_DEBUG("FIRST REL = %p",rp[0]); */

  mtx_output_assign(mtx,rlen,vlen);
  rank = mtx_symbolic_rank(mtx);
  if (rank == 0 ) return 1;     /* nothing to do, eh? */

  /* CONSOLE_DEBUG("FIRST REL = %p",rp[0]); */

  /* lot of whining about dof */
  if (rank < nrow) {
    ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"System is row rank deficient (%d dependent equations)",
            nrow - rank);
  }
  if (rank < ncol) {
    if ( nrow != rank) {
      ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"System is row rank deficient with %d excess columns.",
              ncol - rank);
    } else {
      ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"System has %d degrees of freedom.", ncol - rank);
    }
  }
  if (ncol == nrow) {
    if (ncol != rank) {
      ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"System is (%d) square but rank deficient.",ncol);
    } else {
      ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"System is (%d) square.",ncol);
    }
  }
  if (uppertriangular) {
    mtx_ut_partition(mtx);
  } else {
    mtx_partition(mtx);
  }
  /* copy the block list. there is at least one if rank >=1 as above */
  len = mtx_number_of_blocks(mtx);
  newblocks = ASC_NEW_ARRAY(mtx_region_t,len);
  if (newblocks == NULL) {
    mtx_destroy(mtx);
    return 2;
  }
  for (c = 0 ; c < len; c++) {
    mtx_block(mtx,c,&(newblocks[c]));
  }
  slv_set_solvers_blocks(sys,len,newblocks); /* give it to the system */
  d = slv_get_dofdata(sys);
  d->structural_rank = rank;
  d->n_rows = nrow;
  d->n_cols = ncol;

  /* CONSOLE_DEBUG("FIRST REL = %p",rp[0]); */

  /* the next two lines assume inequalities are un-included. */
#define MINE 1
#if MINE
  d->n_fixed = ncolpfix - ncol;
  d->n_unincluded = nrowpun - nrow;
#else
  d->n_fixed = vlen - ncol;
  d->n_unincluded = rlen - nrow;
#endif
  d->reorder.partition = 1;     /* yes */
  d->reorder.basis_selection = 0;   /* none yet */
  d->reorder.block_reordering = 0;  /* none */

#if SBPDEBUG
  fp = fopen("/tmp/sbp1.plot","w+");
  if (fp !=NULL) {
    mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
    fclose(fp);
  }
#endif

#if MINE
  len = vlenmnv - 1;  /* row of last possible solver variable */
  for (c=ncol; len > c ; ) { /* sort the fixed out of inactive */
    r = mtx_col_to_org(mtx,c);
    if ( var_fixed(vp[r]) && var_active(vp[r])) {
      c++;
    } else {
      mtx_swap_cols(mtx,c,len);
      len--;
    }
  }
#endif

  /* Now, leftover columns are where we want them.
   * That is, unassigned incidence is next to assigned and fixed on right with
   * no columns which do not correspond to variables in the range 0..vlen-1.
   */
  /* there aren't any nonincident vars in the solvers varlist, so
   * we now have the columns in |assigned|dof|fixed|inactive|nonvars| order */

  /* rows 0->rlen-1 should be actual rels, though, and not phantoms.
   * At this point we should have
   * rows - assigned
   *      - unassigned w/incidence
   *      - unincluded and included w/out free incidence and inactive mixed.
   */
  /* sort unincluded relations to bottom of good region */
  /* this is needed because the idiot user may have fixed everything */
  /* in a relation, but we want those rows not to be confused with */
  /* unincluded rows.  we're sort of dropping inequalities on the floor.*/
  len = rlen - 1;  /* row of last possible relation */
  /* sort the inactive out of the unassigned/unincluded */
  for (c=rank; len > c ; ) {
    r = mtx_row_to_org(mtx,c);
    if ( rel_active(rp[r])) {
      c++;
    } else {
      mtx_swap_rows(mtx,c,len);
      len--;
    }
  }

#if MINE
  /* sort the unincluded out of the unassigned */
  len = nrowpun - 1;  /* row of last possible unassigned/unincluded relation */
  for (c=rank; len > c ; ) {
    r = mtx_row_to_org(mtx,c);
    if ( rel_included(rp[r]) ) {
      c++;
    } else {
      mtx_swap_rows(mtx,c,len);
      len--;
    }
  }
#endif


#if SBPDEBUG
  fp = fopen("/tmp/sbp2.plot","w+");
  if (fp !=NULL) {
    mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
    fclose(fp);
  }
#endif
  /* now, at last we have cols jacobian in the order we want the lists to
   * be handed to the solvers. So, we need to reset the var sindex values
   * and reorder the lists.
   */
  if (reindex_vars_from_mtx(sys,0,vlen-1,mtx)) {
    mtx_destroy(mtx);
    return 2;
  }
  /* now, at last we have rows jacobian in the order we want the lists to
   * be handed to the solvers. So, we need to reset the rel sindex values.
   */
  if (reindex_rels_from_mtx(sys,0,rlen-1,mtx)) {
    mtx_destroy(mtx);
    return 2;
  }

  /* CONSOLE_DEBUG("FIRST REL = %p",rp[0]); */

  mtx_destroy(mtx);
  return 0;
}


#if 0 /* code not currently used */
#\ifdef STATIC_HARWELL /* added the backslash so that syntax highlighting behaves */
/* note that you can't statically link to Harwell routines under the terms of the GPL */
extern void mc21b();
extern void mc13emod();
/* returns 0 if ok, OTHERWISE if madness detected.
 * doesn't grok inequalities.
 * CURRENTLY DOESN'T DO WELL WHEN NCOL<NROW
 */
#define SBPDEBUG 0
int slv_block_partition_harwell(slv_system_t sys)
{
#\if SBPDEBUG
  FILE *fp;
#\endif
  struct rel_relation **rp;
  struct rel_relation **rtmp;
  struct rel_relation *rel;
  struct var_variable **vp;
  struct var_variable **vtmp;
  struct var_variable *var;
  const struct var_variable **list;
  mtx_region_t  *newblocks;
  dof_t *d;
  int32 nrow,ncol,size;
  int32 order, rank;
  int32 c,len,vlen,r,v,rlen,col,row;
  int32 licn,rel_tmp_end,rel_count,var_perm_count,var_tmp_end,row_perm_count;
  int32 col_incidence,numnz,row_high,col_high,dummy_ptr,row_ptr,num;
  int32 *dummy_array_ptr;
  var_filter_t vf;
  rel_filter_t rf;
  int32 *icn, *ip, *lenr, *iperm, *iw1, *iw2, *iw3, *arp, *ipnew, *ib;
  int32 *row_perm;

  /* get rel and var info */
  rp = slv_get_solvers_rel_list(sys);
  vp = slv_get_solvers_var_list(sys);
  rlen = slv_get_num_solvers_rels(sys);
  vlen = slv_get_num_solvers_vars(sys);
  if (rlen ==0 || vlen == 0) return 1;
  order = MAX(rlen,vlen);

  /* allocate temp arrays */
  vtmp = (struct var_variable **)ascmalloc(vlen*sizeof(struct var_variable *));
  var = (struct var_variable *)ascmalloc(sizeof(struct var_variable *));
  if (vtmp == NULL || var ==NULL) {
    return 1;
  }
  rtmp = ASC_NEW_ARRAY(struct rel_relation *,rlen);
  rel = ASC_NEW(struct rel_relation);
  if (rtmp == NULL || rel ==NULL) {
    return 1;
  }

  /* set up filters */
  rf.matchbits = (REL_INCLUDED | REL_EQUALITY);
  rf.matchvalue = (REL_INCLUDED | REL_EQUALITY);
  vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED);
  vf.matchvalue = (VAR_INCIDENT | VAR_SVAR);

  /* count rows and cols */
  nrow = slv_count_solvers_rels(sys,&rf);
  ncol = slv_count_solvers_vars(sys,&vf);


  /* count incidence for equality relations that are included
     and active.  Sorts out unincluded and inactive rels to
     the end of the temporary relation list */
  licn = 0;
  rel_tmp_end = rlen -1;
  for (r=0; r < rlen; r++) {
    rel = rp[r];
    if (rel_apply_filter(rel,&rf) && rel_active(rel)) {
      len = rel_n_incidences(rel);
      list = rel_incidence_list(rel);
      for (c=0; c < len; c++) {
    if( var_apply_filter(list[c],&vf) ) {
      licn++;
    }
      }
    } else {
      rtmp[rel_tmp_end] = rp[r];
      rel_tmp_end--;
    }
  }

  /* sort the inactive out of the unincluded and move to end */
  len = rlen -1;
  for(c = rel_tmp_end + 1; len > c; ) {
    rel = rtmp[c];
    if (rel_active(rel)) {
      c++;
    } else {
      rtmp[c] = rtmp[len];
      rtmp[len] = rel;
      len--;
    }
  }

  /* Sort variable list */
  var_perm_count = 0;
  var_tmp_end = vlen - 1;
  for (c = 0; c < vlen; c++) {
    var = vp[c];
    if (var_apply_filter(var,&vf)) {
      vtmp[var_perm_count] = var;
      var_perm_count++;
    } else {
      vtmp[var_tmp_end] = var;
      var_tmp_end--;
    }
  }
  /* sort the inactive out of the unincluded and move to end */
  len = vlen -1;
  for(c = var_tmp_end + 1; len > c; ) {
    var = vtmp[c];
    if (var_active(var) && var_active(var)) {
      c++;
    } else {
      vtmp[c] = vtmp[len];
      vtmp[len] = var;
      len--;
    }
  }

  /* reset solver indicies to make life easier */
  for (c = 0; c < vlen; c++) {
    vp[c] = vtmp[c];
    var_set_sindex(vp[c],c);
  }

  size = MAX(nrow,ncol);
  /* Create vectors for fortran calls */
  icn = ASC_NEW_ARRAY(int32,licn);
  ip = ASC_NEW_ARRAY(int32,size);
  lenr = ASC_NEW_ARRAY(int32,size);
  iperm = ASC_NEW_ARRAY(int32,size);
  iw1 = ASC_NEW_ARRAY(int32,size);
  iw2 = ASC_NEW_ARRAY(int32,size);
  iw3 = ASC_NEW_ARRAY(int32,size);
  arp = ASC_NEW_ARRAY(int32,size);
  ipnew = ASC_NEW_ARRAY(int32,size);
  ib = ASC_NEW_ARRAY(int32,size);
  row_perm = ASC_NEW_ARRAY(int32,size);


/* Fill incidence vectors and place included relations w/o
 * incidence in temporary relation list before the unincluded
 * and inactive relations
 */
  col_incidence = 0;
  row_perm_count = 0;
  for (r=0; r < rlen; r++) {
    rel = rp[r];
    if (rel_apply_filter(rel,&rf) && rel_active(rel)) {
      len = rel_n_incidences(rel);
      if (len > 0) {
    ip[row_perm_count] = col_incidence + 1; /*FORTRAN*/
    lenr[row_perm_count] = len;
    row_perm[row_perm_count] = r;
    row_perm_count++;
    list = rel_incidence_list(rel);
    for (c=0; c < len; c++) {
      if( var_apply_filter(list[c],&vf) ) {
        col = var_sindex(list[c]);
        icn[col_incidence] = col + 1; /*FORTRAN*/
        col_incidence++;
      }
    }
      } else {
    rtmp[rel_tmp_end] = rel;
    rel_tmp_end--;
      }
    }
  }

  licn = col_incidence;
  order = row_perm_count;
  mc21b(&order,icn,&licn,ip,lenr,iperm,&numnz,iw1,arp,iw2,iw3);
  if (order == numnz){
    FPRINTF(stderr,"they are equal\n");
  } else if(order == numnz + 1){
    FPRINTF(stderr,"ADD ONE\n");
  } else {
    FPRINTF(stderr,"no relationship\n");
  }

  if (rank == 0 ) {
    return 1;   /* nothing to do, eh? */
  }
  row_high = nrow - 1;
  col_high = ncol -1;

  /* lot of whining about dof */
  if (rank < nrow) {
    FPRINTF(stderr,"System is row rank deficient (dependent equations)\n");
    row_high = rank - 1;
    /* KHACK: have found that this is executed erroneously */
  }
  if (rank < ncol) {
    if ( nrow != rank) {
      FPRINTF(stderr,"System is row rank deficient with excess columns.\n");
    } else {
      FPRINTF(stderr,"System has degrees of freedom.\n");
    }
  }
  if (ncol == nrow) {
    FPRINTF(stderr,"System is (%d) square",ncol);
    if (ncol != rank) {
      FPRINTF(stderr,"but rank deficient.\n");
    } else {
      FPRINTF(stderr,".\n");
    }
  }

  dummy_array_ptr = arp;
  arp = lenr;
  lenr = dummy_array_ptr;

  /* Fix row pointers for assigned relations. Put unassigned
     onto temporary relation list (before included w/o incidence)
     and put dummy dense rows into ipnew */
  dummy_ptr = licn + 1; /* FORTRAN */
  for (row = 0; row < order; row++) {
    row_ptr = iperm[row] - 1; /* FORTRAN */
    if (row_ptr != -1) { /* indicates unassigned row */
      ipnew[row_ptr] = ip[row];
      lenr[row_ptr] = arp[row];
    } else {
      ipnew[row_ptr] = dummy_ptr;
      dummy_ptr += ncol;
      lenr[row_ptr] = ncol;
      /* unassigned before fixed*/
      rtmp[rel_tmp_end] = rp[row_perm[row_ptr]];
      rel_tmp_end--;
    }
  }
  for (row = order; row < ncol; row++) {
    ipnew[row] = dummy_ptr;
    dummy_ptr += ncol;
    lenr[row] = ncol;
  }

  mc13emod(&size,icn,&licn,ipnew,lenr,arp,ib,&num,iw1,iw2,iw3);


  /* copy the block list. there is at least one if rank >=1 as above */
  len = num;
  newblocks = ASC_NEW_ARRAY(mtx_region_t,len);
  if (newblocks == NULL) {
    return 2;
  }
  /* set up locations of block corners with the row and column
     orderings which will be set soon */
  for (c = 0 ; c < len; c++) {
    newblocks[c].row.low = ib[c]-1;
    newblocks[c].col.low = newblocks[c].row.low;
  }
  for (c = 0 ; c < len -1; c++) {
    newblocks[c].row.high = newblocks[c+1].row.low -1;
    newblocks[c].col.high = newblocks[c].row.high;
  }
  newblocks[len - 1].row.high = row_high;
  newblocks[len - 1].col.high = col_high;

  slv_set_solvers_blocks(sys,len,newblocks); /* give it to the system */
  d = slv_get_dofdata(sys);
  d->structural_rank = rank;
  d->n_rows = nrow;
  d->n_cols = ncol;
  /* the next two lines assume only free and fixed incident solvervar
   * appear in the solvers var list and inequalities are unincluded.
   */
  d->n_fixed = vlen - ncol;
  d->n_unincluded = rlen - nrow;
  d->reorder.partition = 1;     /* yes */
  d->reorder.basis_selection = 0;   /* none yet */
  d->reorder.block_reordering = 0;  /* none */

  for (c = 0; c < ncol; c++) {
    v = arp[c] -1; /* FORTRAN */
    vtmp[c] = vp[v];
  }
  for (c = 0; c < vlen; c++) {
    vp[c] = vtmp[c];
    var_set_sindex(vp[c],c);
  }

  rel_count = 0;
  for (c = 0; c < size; c++) {
    r = arp[c] - 1;
    if (r < order){ /* weed out fake rows */
      r = iperm[r] - 1;
      if (r != -1) { /* unassigned already in rtmp */
    rtmp[rel_count] = rp[row_perm[r]];
    rel_count++;
      }
    }
  }
  for (c = 0; c < rlen; c++) {
    rp[c] = rtmp[c];
    rel_set_sindex(rp[c],c);
  }
  ascfree(vtmp);
  ascfree(rtmp);
  ascfree(icn);
  ascfree(ip);
  ascfree(lenr);
  ascfree(iperm);
  ascfree(iw1);
  ascfree(iw2);
  ascfree(iw3);
  ascfree(arp);
  ascfree(ipnew);
  ascfree(ib);
  return 0;
}
#\endif /*STATIC_HARWELL*/
#endif /* 0 */


void slv_sort_rels_and_vars(slv_system_t sys,
                int32 *rel_count,
                int32 *var_count)
{
  struct rel_relation **rp;
  struct rel_relation **rtmp;
  struct rel_relation *rel;
  struct var_variable **vp;
  struct var_variable **vtmp;
  struct var_variable *var;
  int32 nrow,ncol,rlen,vlen,order,rel_tmp_end;
  int32 r,c,var_tmp_end,len;
  var_filter_t vf;
  rel_filter_t rf;

  /* get rel and var info */
  rp = slv_get_solvers_rel_list(sys);
  vp = slv_get_solvers_var_list(sys);
  rlen = slv_get_num_solvers_rels(sys);
  vlen = slv_get_num_solvers_vars(sys);
  if (rlen ==0 || vlen == 0) return;
  order = MAX(rlen,vlen);

  assert(var_count != NULL && rel_count != NULL);
  *var_count = *rel_count = -1;

  /* allocate temp arrays */
  vtmp = (struct var_variable **)ascmalloc(vlen*sizeof(struct var_variable *));
  if (vtmp == NULL) {
    return;
  }
  rtmp = ASC_NEW_ARRAY(struct rel_relation *,rlen);
  if (rtmp == NULL) {
    ascfree(vtmp);
    return;
  }

  /* set up filters */
  rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
  rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
  vf.matchbits = (VAR_INCIDENT | VAR_SVAR | VAR_FIXED | VAR_ACTIVE);
  vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);

  /* count rows and cols */
  nrow = slv_count_solvers_rels(sys,&rf);
  ncol = slv_count_solvers_vars(sys,&vf);


  /* Sort out unincluded and inactive rels to
   * the end of the temporary relation list.
   */
  *rel_count = 0;
  rel_tmp_end = rlen -1;
  for (r=0; r < rlen; r++) {
    rel = rp[r];
    if (rel_apply_filter(rel,&rf)) {
      rtmp[*rel_count] = rel;
      (*rel_count)++;
    } else {
      rtmp[rel_tmp_end] = rel;
      rel_tmp_end--;
    }
  }

  /* sort the inactive out of the unincluded and move to end */
  len = rlen -1;
  for(c = rel_tmp_end + 1; len > c; ) {
    rel = rtmp[c];
    if (rel_active(rel)) {
      c++;
    } else {
      rtmp[c] = rtmp[len];
      rtmp[len] = rel;
      len--;
    }
  }

  /* Sort variable list */
  *var_count = 0;
  var_tmp_end = vlen - 1;
  for (c = 0; c < vlen; c++) {
    var = vp[c];
    if (var_apply_filter(var,&vf)) {
      vtmp[*var_count] = var;
      (*var_count)++;
    } else {
      vtmp[var_tmp_end] = var;
      var_tmp_end--;
    }
  }
  /* sort the inactive out of the unincluded and move to end */
  len = vlen -1;
  for(c = var_tmp_end + 1; len > c; ) {
    var = vtmp[c];
    if (var_active(var) && var_active(var)) {
      c++;
    } else {
      vtmp[c] = vtmp[len];
      vtmp[len] = var;
      len--;
    }
  }

  /* reset solver indicies */
  for (c = 0; c < vlen; c++) {
    vp[c] = vtmp[c];
    var_set_sindex(vp[c],c);
  }

  for (c = 0; c < rlen; c++) {
    rp[c] = rtmp[c];
    rel_set_sindex(rp[c],c);
  }
  ascfree(vtmp);
  ascfree(rtmp);
  return;
}


int slv_block_unify(slv_system_t sys)
{
  dof_t *d;
  const mtx_block_t *mbt;
  mtx_region_t *newblocks;

  if (sys==NULL) return 1;

  d = slv_get_dofdata(sys);
  mbt = slv_get_solvers_blocks(sys);
  assert(d!=NULL && mbt!=NULL);
  if (d->structural_rank && mbt->nblocks >1) {
    newblocks = ASC_NEW(mtx_region_t);
    if (newblocks == NULL) return 2;
    newblocks->row.low = newblocks->col.low = 0;
    newblocks->row.high = d->structural_rank - 1;
    newblocks->col.high = d->n_cols - 1;
    if ( newblocks->col.high < 0 || newblocks->row.high < 0) {
      ascfree(newblocks);
      return 1;
    }
    slv_set_solvers_blocks(sys,1,newblocks);
  }
  return 0;
}

int slv_set_up_block(slv_system_t sys,int bnum)
{
  struct rel_relation **rp;
  struct var_variable **vp;
  mtx_region_t reg;
  const mtx_block_t *b;
  int32 c,vlen,rlen;
  var_filter_t vf;
  rel_filter_t rf;

  if (sys==NULL) return 1;
  rlen = slv_get_num_solvers_rels(sys);
  vlen = slv_get_num_solvers_vars(sys);
  if (rlen ==0 || vlen == 0) return 1;

  rp = slv_get_solvers_rel_list(sys);
  vp = slv_get_solvers_var_list(sys);
  assert(rp!=NULL);
  assert(vp!=NULL);

  rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
  rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
  vf.matchbits =(VAR_INCIDENT |VAR_SVAR | VAR_FIXED |VAR_INBLOCK | VAR_ACTIVE);
  vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_INBLOCK | VAR_ACTIVE);

  b = slv_get_solvers_blocks(sys);
  assert(b!=NULL);
  if (bnum <0 || bnum >= b->nblocks || b->block == NULL) return 1;
  reg = b->block[bnum];
  for (c=reg.col.low; c<=reg.col.high; c++) {
    var_set_in_block(vp[c],1);
  }
  for (c=reg.row.low; c<=reg.row.high; c++) {
    rel_set_in_block(rp[c],1);
  }
  return 0;
}

/* returns 0 if ok, OTHERWISE if madness detected.
 */
int slv_spk1_reorder_block(slv_system_t sys,int bnum,int transpose)
{
  struct rel_relation **rp;
  struct var_variable **vp;
  mtx_region_t reg;
  const mtx_block_t *b;
  mtx_matrix_t mtx;
  mtx_coord_t coord;
  int32 c,vlen,rlen;
  var_filter_t vf;
  rel_filter_t rf;
  dof_t *d;

  if (sys==NULL) return 1;
  rlen = slv_get_num_solvers_rels(sys);
  vlen = slv_get_num_solvers_vars(sys);
  if (rlen ==0 || vlen == 0) return 1;

  rp = slv_get_solvers_rel_list(sys);
  vp = slv_get_solvers_var_list(sys);
  assert(rp!=NULL);
  assert(vp!=NULL);
  rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
  rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
  vf.matchbits =(VAR_INCIDENT |VAR_SVAR | VAR_FIXED |VAR_INBLOCK | VAR_ACTIVE);
  vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_INBLOCK | VAR_ACTIVE);

  mtx = mtx_create();
  mtx_set_order(mtx,MAX(rlen,vlen));
  b = slv_get_solvers_blocks(sys);
  assert(b!=NULL);
  if (bnum <0 || bnum >= b->nblocks || b->block == NULL) return 1;
  reg = b->block[bnum];
  for (c=reg.col.low; c<=reg.col.high; c++) {
    var_set_in_block(vp[c],1);
  }
  for (c=reg.row.low; c<=reg.row.high; c++) {
    rel_set_in_block(rp[c],1);
  }
  if (reg.row.low != reg.col.low || reg.row.high != reg.col.high) {
    return 1; /* must be square */
    /* could also enforce minsize 3x3, but someone my call with a
     * partitionable region, so don't want to.
     */
  }

  if (slv_std_make_incidence_mtx(sys,mtx,&vf,&rf)) {
    FPRINTF(stderr,
      "slv_spk1_reorder_block: failure in creating incidence matrix.\n");
    mtx_destroy(mtx);
    return 1;
  }
  /* verify that block has no empty columns, though not checking diagonal */
  for (c = reg.row.low; c <= reg.row.high; c++) {
    coord.col = mtx_FIRST;
    coord.row = c;
    if (mtx_next_in_row(mtx,&coord,mtx_ALL_COLS), coord.col == mtx_LAST) {
      mtx_destroy(mtx);
      FPRINTF(stderr, "slv_spk1_reorder_block: empty row (%d) found.\n",c);
      return 1;
    }
    coord.row = mtx_FIRST;
    coord.col = c;
    if (mtx_next_in_col(mtx,&coord,mtx_ALL_ROWS), coord.row == mtx_LAST) {
      FPRINTF(stderr, "slv_spk1_reorder_block: empty col (%d) found.\n",c);
      mtx_destroy(mtx);
      return 1;
    }
  }
  if (transpose) {
    mtx_reorder(mtx,&reg,mtx_TSPK1);
  } else {
    mtx_reorder(mtx,&reg,mtx_SPK1);
  }
  if (reindex_vars_from_mtx(sys,reg.col.low,reg.col.high,mtx)) {
    mtx_destroy(mtx);
    return 2;
  }
  if (reindex_rels_from_mtx(sys,reg.row.low,reg.row.high,mtx)) {
    mtx_destroy(mtx);
    return 2;
  }
  d = slv_get_dofdata(sys);
  d->reorder.block_reordering = 1;  /* spk1 */

  mtx_destroy(mtx);
  return 0;
}

/*
    A function to be called on the leftover diagonal blocks of mr_bisect.
    This function should be thoroughly investigated, which it has not been.
*/
static int tear_subreorder(slv_system_t server,
                           mtx_matrix_t mtx, mtx_region_t *reg)
{
  assert(mtx != NULL && reg !=NULL);
  (void)server;
  if (g_blockmethod==mtx_UNKNOWN) {
    mtx_reorder(mtx,reg,mtx_SPK1);
  } else {
    mtx_reorder(mtx,reg,g_blockmethod);
  }
  return 0;
}

int slv_tear_drop_reorder_block(slv_system_t sys, int32 bnum,
        int32 cutoff, int two,
        enum mtx_reorder_method blockmethod
){
  struct rel_relation **rp;
  struct var_variable **vp;
  const mtx_block_t *b;
  dof_t *d;
  mr_reorder_t *mrsys;
  mtx_matrix_t mtx;
  mtx_region_t reg;
  mtx_coord_t coord;
  int32 c,vlen,rlen,modcount;
  var_filter_t vf;
  rel_filter_t rf;

  if (sys==NULL) return 1;
  rlen = slv_get_num_solvers_rels(sys);
  vlen = slv_get_num_solvers_vars(sys);
  if (rlen ==0 || vlen == 0) return 1;

  rp = slv_get_solvers_rel_list(sys);
  vp = slv_get_solvers_var_list(sys);
  assert(rp!=NULL);
  assert(vp!=NULL);
  rf.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK | REL_ACTIVE);
  rf.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_INBLOCK| REL_ACTIVE);
  vf.matchbits =(VAR_INCIDENT |VAR_SVAR | VAR_FIXED |VAR_INBLOCK | VAR_ACTIVE);
  vf.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_INBLOCK | VAR_ACTIVE);

  mtx = mtx_create();
  mtx_set_order(mtx,MAX(rlen,vlen));
  b = slv_get_solvers_blocks(sys);
  assert(b!=NULL); /* probably shouldn't be an assert here ... */
  if (bnum <0 || bnum >= b->nblocks || b->block == NULL) {
    mtx_destroy(mtx);
    return 1;
  }
  reg = b->block[bnum];
  /* do the flag setting even if we don't do the reorder */
  for (c=reg.col.low; c<=reg.col.high; c++) {
    var_set_in_block(vp[c],1);
  }
  for (c=reg.row.low; c<=reg.row.high; c++) {
    rel_set_in_block(rp[c],1);
  }
  if (reg.row.low != reg.col.low || reg.row.high != reg.col.high) {
    mtx_destroy(mtx);
    return 1; /* must be square */
  }
  if (reg.row.high - reg.row.low < 3) {
    mtx_destroy(mtx);
    return 0; /* must be 3x3 or bigger to have any effect */
  }

  if (slv_std_make_incidence_mtx(sys,mtx,&vf,&rf)) {
    FPRINTF(stderr,
      "slv_tear_drop_reorder_block: failure in creating incidence matrix.\n");
    mtx_destroy(mtx);
    return 1;
  }
  /* verify that block has no empty columns, though not checking diagonal */
  for (c = reg.row.low; c <= reg.row.high; c++) {
    coord.col = mtx_FIRST;
    coord.row = c;
    if (mtx_next_in_row(mtx,&coord,mtx_ALL_COLS), coord.col == mtx_LAST) {
      mtx_destroy(mtx);
      FPRINTF(stderr, "slv_tear_drop_reorder_block: empty row (%d) found.\n",c);
      return 1;
    }
    coord.row = mtx_FIRST;
    coord.col = c;
    if (mtx_next_in_col(mtx,&coord,mtx_ALL_ROWS), coord.row == mtx_LAST) {
      FPRINTF(stderr, "slv_tear_drop_reorder_block: empty col (%d) found.\n",c);
      mtx_destroy(mtx);
      return 1;
    }
  }

  modcount = slv_get_num_models(sys);
  mrsys = mr_reorder_create(sys,mtx,modcount);
  if (mrsys==NULL) return 1;
  mrsys->cutoff = cutoff;
  g_blockmethod = blockmethod;
  if (!two) {
    mr_bisect_partition(mrsys,&reg,1,tear_subreorder);
  } else {
    mr_bisect_partition2(mrsys,&reg,1,tear_subreorder);
  }
#if 1
  FPRINTF(stderr,"Model-based reordering: (tear-style = %d)\n",two);
  FPRINTF(stderr,"Min. tearable size:\t%d\n",mrsys->cutoff);
  FPRINTF(stderr,"Tears:\t\t%d\n",mrsys->ntears);
  FPRINTF(stderr,"Partitionings:\t%d\n",mrsys->nblts);
#endif
  reg = b->block[bnum]; /* bisect likely munged reg */
  if (reindex_vars_from_mtx(sys,reg.col.low,reg.col.high,mtx)) {
    mtx_destroy(mtx);
    mr_reorder_destroy(mrsys);
    return 2;
  }
  if (reindex_rels_from_mtx(sys,reg.row.low,reg.row.high,mtx)) {
    mtx_destroy(mtx);
    mr_reorder_destroy(mrsys);
    return 2;
  }
  d = slv_get_dofdata(sys);
  d->reorder.block_reordering = 2;  /* tear_drop_baa */

  mtx_destroy(mtx);
  mr_reorder_destroy(mrsys);
  return 0;
}

int slv_insure_bounds(slv_system_t sys,int32 lo,int32 hi, FILE *mif)
{
  real64 val,low,high;
  int32 c,nchange=0;
  struct var_variable *var, **vp;

  vp = slv_get_solvers_var_list(sys);
  if (vp==NULL) return -1;
  for (c= lo; c <= hi; c++) {
    var = vp[c];
    low = var_lower_bound(var);
    high = var_upper_bound(var);
    val = var_value(var);
    if( low > high ) {
      if (mif!=NULL) {
        FPRINTF(mif,"Bounds for variable ");
        var_write_name(sys,var,mif);
        FPRINTF(mif," are inconsistent [%g,%g].\n",low,high);
        FPRINTF(mif,"Bounds will be swapped.\n");
      }
      var_set_upper_bound(var, low);
      var_set_lower_bound(var, high);
      low = var_lower_bound(var);
      high = var_upper_bound(var);
      nchange++;
    }

    if( low > val ) {
      if (mif!=NULL) {
        FPRINTF(mif,"Variable ");
        var_write_name(sys,var,mif);
        FPRINTF(mif," was initialized below its lower bound.\n");
        FPRINTF(mif,"It will be moved to its lower bound.\n");
      }
      var_set_value(var, low);
    } else {
      if( val > high ) {
        if (mif!=NULL) {
          FPRINTF(mif,"Variable ");
          var_write_name(sys,var,mif);
          FPRINTF(mif," was initialized above its upper bound.\n");
          FPRINTF(mif,"It will be moved to its upper bound.\n");
        }
        var_set_value(var, high);
        nchange++;
      }
    }
  }
  return nchange;
}


void slv_check_bounds(const slv_system_t sys,int32 lo,int32 hi,
                      FILE *mif,const char *label)
{
  real64 val,low,high;
  int32 c,len;
  struct var_variable *var, **vp;
  static char defaultlabel[] = "";

  if (label==NULL) label = defaultlabel;
  if (sys==NULL || mif == NULL) return;
  vp = slv_get_solvers_var_list(sys);
  if (vp==NULL) return;
  len = slv_get_num_solvers_vars(sys);
  if (lo > len || hi > len) {
    FPRINTF(stderr,"slv_check_bounds miscalled\n");
    return;
  }
  for (c= lo; c <= hi; c++) {
    var = vp[c];
    low = var_lower_bound(var);
    high = var_upper_bound(var);
    val = var_value(var);
    if( low > high ) {
      FPRINTF(mif,"Bounds for %s variable ",label);
      var_write_name(sys,var,mif);
      FPRINTF(mif,"\nare inconsistent [%g,%g].\n",low,high);
    }

    if( low > val ) {
      FPRINTF(mif,"%s variable ",label);
      var_write_name(sys,var,mif);
      FPRINTF(mif,"\nwas initialized below its lower bound.\n");
    } else {
      if( val > high ) {
        FPRINTF(mif,"%s variable ",label);
        var_write_name(sys,var,mif);
        FPRINTF(mif," was initialized above its upper bound.\n");
      }
    }
  }
  return;
}

/*------------------------------------------------------------------------------
  SOLVER REGISTRATION
*/

/* rewrote this stuff to get rid of all the #ifdefs -- JP */

struct StaticSolverRegistration{
    int is_active;
    const char *name;
    SlvRegistration *regfunc;
};

/*
    The names here are only used to provide information in the case where
    solver registration fails. The definitive solver names are in the slv*.c
    files.
*/
static const struct StaticSolverRegistration slv_reg[]={
/*  {0,"SLV",&slv0_register} */
/*  ,{0,"MINOS",&slv1_register} */
    {HAVE_QRSLV,"QRSLV",&slv3_register}
/*  ,{0,"CSLV",&slv4_register} */
/*  ,{0,"LSSLV",&slv5_register} */
/*  ,{0,"MPS",&slv6_register} */
/*  ,{0,"NGSLV",&slv7_register} */
/*  ,{0,"OPTSQP",&slv2_register} */
    ,{HAVE_CONOPT,"CONOPT",&slv8_register}
    ,{HAVE_CMSLV,"CMSLV",&slv9_register}
/*  ,{0,"LRSLV",&slv9a_register} */
    ,{0,NULL,NULL}
};
    
int SlvRegisterStandardClients(void){
    int nclients = 0;
    int newclient=0;
    int error;
    int i;
    /* CONSOLE_DEBUG("REGISTERING STANDARD SOLVER ENGINES"); */
    for(i=0;slv_reg[i].name!=NULL;++i){
        if(slv_reg[i].is_active && slv_reg[i].regfunc){
            error = slv_register_client(slv_reg[i].regfunc,NULL,NULL,&newclient);
            if(error){
                ERROR_REPORTER_HERE(ASC_PROG_ERR
                    ,"Unable to register solver '%s' (error %d)."
                    ,slv_reg[i].name,error
                );
            }else{
                CONSOLE_DEBUG("Solver '%s' registered OK",slv_solver_name(newclient));
                nclients++;
            }
        }else{
            CONSOLE_DEBUG("Solver '%s' was not compiled.",slv_reg[i].name);
        }
    }
  return nclients;
}

/*------------------------------------------------------------------------------
  OUTPUT ASSIGNMENT AND PARTITIONG IN LOGICAL RELATIONS
*/

#define MLIMDEBUG 0
/*
 * Populates a matrix according to the sys solvers_dvars, solvers_logrels
 * lists and the filters given. mtx given must be created (not null)
 * and with order >= MAX( slv_get_num_solvers_logrels(sys),
 *                        slv_get_num_solvers_dvars(sys));
 * returns 0 if went ok.
 */
int slv_std_make_log_incidence_mtx(slv_system_t sys, mtx_matrix_t mtx,
                                   dis_filter_t *dvf,logrel_filter_t *lrf)
{
#if MLIMDEBUG
  FILE *fp;
#endif
  int32 lr, lrlen,ord;
  struct logrel_relation **lrp;

  if (sys==NULL || mtx == NULL || dvf == NULL || lrf == NULL) {
    FPRINTF(stderr,"make_log_incidence called with null\n");
    return 1;
  }
  lrp = slv_get_solvers_logrel_list(sys);
  assert(lrp!=NULL);
  lrlen = slv_get_num_solvers_logrels(sys);
  ord = MAX(lrlen,slv_get_num_solvers_dvars(sys));
  if (ord > mtx_order(mtx)) {
    FPRINTF(stderr,"make_incidence called with undersized matrix\n");
    return 2;
  }
  for (lr=0; lr < lrlen; lr++) {
    if (logrel_apply_filter(lrp[lr],lrf)) {
      logrelman_get_incidence(lrp[lr],dvf,mtx);
    }
  }
#if MLIMDEBUG
  fp = fopen("/tmp/mim.plot","w+");
  if (fp !=NULL) {
    mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
    fclose(fp);
  }
#endif
  return 0;
}

#if USEDCODE

/* returns 0 if successful, 1 if insufficient memory. Does not change
 * the data in mtx. Orders the solvers_dvar list of the system to
 * match the permutation on the given mtx. It is assumed that
 * the org cols of mtx == dvar list position. Only the
 * dvars in range lo to hi of the dvar list are affected and
 * only these dvars should appear in the same org column range of
 * the mtx. This should not be called on blocks less than 3x3.
 */
static int reindex_dvars_from_mtx(slv_system_t sys, int32 lo, int32 hi,
                                  const mtx_matrix_t mtx)
{
  struct dis_discrete **dvtmp, **dvp;
  int32 c,v,vlen;

  if (lo >= hi +1) return 0; /* job too small */
  dvp = slv_get_solvers_dvar_list(sys);
  vlen = slv_get_num_solvers_dvars(sys);
  /* on dvtmp we DONT have the terminating null */
  dvtmp = (struct dis_discrete **)
                                ascmalloc(vlen*sizeof(struct dis_discrete *));
  if (dvtmp == NULL) {
    return 1;
  }
  /* copy pointers to dvtmp in order desired and copy back rather than sort */
  for (c=lo;c<=hi;c++) {
    v = mtx_col_to_org(mtx,c);
    dvtmp[c] = dvp[v];
  }
  /* copying back and re-sindexing */
  for (c=lo;c<=hi;c++) {
    dvp[c] = dvtmp[c];
    dis_set_sindex(dvp[c],c);
  }
  ascfree(dvtmp);
  return 0;
}

/* returns 0 if successful, 1 if insufficient memory. Does not change
 * the data in mtx. Orders the solvers_logrel list of the system to
 * match the permutation on the given mtx. It is assumed that
 * the org rows of mtx == logrel list position. Only the
 *logrels in range lo to hi of the logrel list are affected and
 * only these logrels should appear in the same org row range of
 * the input mtx. This should not be called on blocks less than 3x3.
 */
static int reindex_logrels_from_mtx(slv_system_t sys, int32 lo, int32 hi,
                                    const mtx_matrix_t mtx)
{
  struct logrel_relation **lrtmp, **lrp;
  int32 c,v,rlen;

  lrp = slv_get_solvers_logrel_list(sys);
  rlen = slv_get_num_solvers_logrels(sys);
  /* on lrtmp we DONT have the terminating null */
  lrtmp = (struct logrel_relation **)
                             ascmalloc(rlen*sizeof(struct logrel_relation *));
  if (lrtmp == NULL) {
    return 1;
  }
  /* copy pointers to lrtmp in order desired and copy back rather than sort.
   * do this only in the row range of interest. */
  for (c=lo;c<=hi;c++) {
    v = mtx_row_to_org(mtx,c);
#if RIDEBUG
    FPRINTF(stderr,"Old logrel sindex(org) %d becoming sindex(cur) %d\n",v,c);
#endif
    lrtmp[c] = lrp[v];
  }
  /* copying back and re-sindexing */
  for (c=lo;c<=hi;c++) {
    lrp[c] = lrtmp[c];
    logrel_set_sindex(lrp[c],c);
  }
  ascfree(lrtmp);
  return 0;
}

#endif /* USEDCODE */

/*
 * returns 0 if ok, OTHERWISE if madness detected.
 */
int slv_log_block_partition(slv_system_t sys)
{
#if SLBPDEBUG
  FILE *fp;
#endif
  struct logrel_relation **lrp;
  struct dis_discrete **dvp;
  mtx_region_t  *newblocks;
  dof_t *d;
  int32 nrow,ncol;
  int32 ncolpfix, nrowpun, vlenmnv;
  mtx_matrix_t mtx;
  int32 order, rank;
  int32 c,len,dvlen,r,lrlen;
  dis_filter_t dvf;
  logrel_filter_t lrf;

  lrp = slv_get_solvers_logrel_list(sys);
  dvp = slv_get_solvers_dvar_list(sys);
  lrlen = slv_get_num_solvers_logrels(sys);
  dvlen = slv_get_num_solvers_dvars(sys);
  if (lrlen ==0 || dvlen == 0) return 1;
  order = MAX(lrlen,dvlen);

  lrf.matchbits = (LOGREL_ACTIVE);
  lrf.matchvalue = (LOGREL_ACTIVE);
  dvf.matchbits = (DIS_INCIDENT | DIS_BVAR  | DIS_ACTIVE);
  dvf.matchvalue = (DIS_INCIDENT | DIS_BVAR  | DIS_ACTIVE);

  /* nrow plus unincluded logrels */
  nrowpun = slv_count_solvers_logrels(sys,&lrf);
  /* ncol plus fixed discrete vars */
  ncolpfix = slv_count_solvers_dvars(sys,&dvf);

  dvf.matchbits = (DIS_BVAR);
  dvf.matchvalue = 0;

  /* dvlen minus non boolean vars */
  vlenmnv = dvlen - slv_count_solvers_dvars(sys,&dvf);

  lrf.matchbits = (LOGREL_INCLUDED | LOGREL_ACTIVE);
  lrf.matchvalue = (LOGREL_INCLUDED | LOGREL_ACTIVE);
  dvf.matchbits = (DIS_INCIDENT | DIS_BVAR | DIS_FIXED | DIS_ACTIVE);
  dvf.matchvalue = (DIS_INCIDENT | DIS_BVAR | DIS_ACTIVE);

  nrow = slv_count_solvers_logrels(sys,&lrf);
  ncol = slv_count_solvers_dvars(sys,&dvf);

  mtx = slv_get_sys_mtx(sys);
  mtx_set_order(mtx,order);

  if (slv_std_make_log_incidence_mtx(sys,mtx,&dvf,&lrf)) {
    FPRINTF(stderr,
      "slv_log_block_partition: failure in creating incidence matrix.\n");
    mtx_destroy(mtx);
    return 1;
  }

  mtx_output_assign(mtx,lrlen,dvlen);
  rank = mtx_symbolic_rank(mtx);
  if (rank == 0 ) return 1;
  if (rank < nrow) {
    FPRINTF(stderr,"rank<nrow in slv_log_block_partition\n");
    FPRINTF(stderr,"dependent logical relations ? \n");
  }
  if (rank < ncol) {
    if ( nrow != rank) {
      FPRINTF(stderr,"rank<ncol and nrow!=rank in slv_log_block_partition.\n");
      FPRINTF(stderr,"Excess of columns ? \n");
    } else {
      FPRINTF(stderr,"rank<ncol but nrow==rank in slv_log_block_partition.\n");
      FPRINTF(stderr,"Degrees of freedom ? \n");
    }
  }
  if (ncol == nrow) {
    if (ncol == rank) {
      FPRINTF(stderr,"\n");
      FPRINTF(stderr,
             "System of logical relations does not need Logic Inference \n");
    }
    if (ncol != rank) {
      FPRINTF(stderr,"but ncol!=rank.\n");
      FPRINTF(stderr,"Rank deficient ? \n");
    } else {
      FPRINTF(stderr,"\n");
    }
  }
  mtx_partition(mtx);
  /* copy the block list. there is at least one if rank >=1 as above */
  len = mtx_number_of_blocks(mtx);
  newblocks = ASC_NEW_ARRAY(mtx_region_t,len);
  if (newblocks == NULL) {
    mtx_destroy(mtx);
    return 2;
  }
  for (c = 0 ; c < len; c++) {
    mtx_block(mtx,c,&(newblocks[c]));
  }
  slv_set_solvers_log_blocks(sys,len,newblocks); /* give it to the system */
  d = slv_get_log_dofdata(sys);
  d->structural_rank = rank;
  d->n_rows = nrow;
  d->n_cols = ncol;
  d->n_fixed = ncolpfix - ncol;
  d->n_unincluded = nrowpun - nrow;
  d->reorder.partition = 1;     /* yes */
  d->reorder.basis_selection = 0;   /* none yet */
  d->reorder.block_reordering = 0;  /* none */

#if SLBPDEBUG
  fp = fopen("/tmp/sbp1.plot","w+");
  if (fp !=NULL) {
    mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
    fclose(fp);
  }
#endif

  len = vlenmnv - 1;  /* row of last possible boolean variable */
  for (c=ncol; len > c ; ) { /* sort the fixed out of inactive */
    r = mtx_col_to_org(mtx,c);
    if ( dis_fixed(dvp[r]) && dis_active(dvp[r])) {
      c++;
    } else {
      mtx_swap_cols(mtx,c,len);
      len--;
    }
  }

  /* Now, leftover columns are where we want them.
   * That is, unassigned incidence is next to assigned and fixed on right with
   * no columns which do not correspond to dis vars in the range 0..dvlen-1.
   */
  /* there aren't any nonincident dvars in the solvers dvarlist, so
   * we now have the columns in |assigned|dof|fixed|inactive|nonvars| order */

  /* rows 0->lrlen-1 should be actual logrels, though, and not phantoms.
   * At this point we should have
   * rows - assigned
   *      - unassigned w/incidence
   *      - unincluded and included w/out free incidence and inactive mixed.
   */

  /* sort unincluded logrelations to bottom of good region */
  /* this is needed because the idiot user may have fixed everything */
  /* in a logrelation, but we want those rows not to be confused with */
  /* unincluded rows.*/
  len = lrlen - 1;  /* row of last possible relation */
  /* sort the inactive out of the unassigned/unincluded */
  for (c=rank; len > c ; ) {
    r = mtx_row_to_org(mtx,c);
    if ( logrel_active(lrp[r])) {
      c++;
    } else {
      mtx_swap_rows(mtx,c,len);
      len--;
    }
  }

  /* sort the unincluded out of the unassigned */
  len = nrowpun - 1;  /* row of last possible unassigned/unincluded relation */
  for (c=rank; len > c ; ) {
    r = mtx_row_to_org(mtx,c);
    if ( logrel_included(lrp[r]) ) {
      c++;
    } else {
      mtx_swap_rows(mtx,c,len);
      len--;
    }
  }

#if SLBPDEBUG
  fp = fopen("/tmp/sbp2.plot","w+");
  if(fp !=NULL){
    mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
    fclose(fp);
  }
#endif

  /* These functions are called for slv_block_partitioning, but I do not
   * know if I need them here:
   *
   * now, at last we have cols in the order we want the lists to
   * be handed to the solver. So, we need to reset the dvar sindex values
   * and reorder the lists.
   */

#if USEDCODE
   if (reindex_dvars_from_mtx(sys,0,dvlen-1,mtx)) {
     mtx_destroy(mtx);
     return 2;
   }
#endif /* USEDCODE */

  /*
   * now, at last we have rows jacobian in the order we want the lists to
   * be handed to the solvers. So, we need to reset the rel sindex values.
   */

#if USEDCODE
   if (reindex_logrels_from_mtx(sys,0,lrlen-1,mtx)) {
     mtx_destroy(mtx);
     return 2;
   }
#endif /* USEDCODE */

  return 0;
}
