cddlp.c~

CheckMate is a MATLAB-based tool for modeling, simulating and investigating properties of hybrid dyn

void dd_WriteTableau(FILE *f,dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
  dd_colindex nbindex,dd_rowindex bflag)
/* Write the tableau  A.T   */
{
  dd_colrange j;
  dd_rowrange i;
  mytype x;
  
  dd_init(x);
  fprintf(f," %ld  %ld  real\n",m_size,d_size);
  fprintf(f,"          |");
  for (j=1; j<= d_size; j++) {
    fprintf(f," %ld",nbindex[j]);
  } fprintf(f,"\n");
  for (j=1; j<= d_size+1; j++) {
    fprintf(f," ----");
  } fprintf(f,"\n");
  for (i=1; i<= m_size; i++) {
    fprintf(f," %3ld(%3ld) |",i,bflag[i]);  
    for (j=1; j<= d_size; j++) {
      dd_TableauEntry(&x,m_size,d_size,A,T,i,j);
      dd_WriteNumber(f,x);
    }
    fprintf(f,"\n");
  }
  fprintf(f,"end\n");
  dd_clear(x);
}

void dd_WriteSignTableau(FILE *f,dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
  dd_colindex nbindex,dd_rowindex bflag)
/* Write the sign tableau  A.T   */
{
  dd_colrange j;
  dd_rowrange i;
  mytype x;
  
  dd_init(x);
  fprintf(f," %ld  %ld  real\n",m_size,d_size);
  fprintf(f,"          |");
  for (j=1; j<= d_size; j++) {
    fprintf(f,"%3ld",nbindex[j]);
  } fprintf(f,"\n  ------- | ");
  for (j=1; j<= d_size; j++) {
    fprintf(f,"---");
  } fprintf(f,"\n");
  for (i=1; i<= m_size; i++) {
    fprintf(f," %3ld(%3ld) |",i,bflag[i]);
    for (j=1; j<= d_size; j++) {
      dd_TableauEntry(&x,m_size,d_size,A,T,i,j);
      if (dd_Positive(x)) fprintf(f, "  +");
      else if (dd_Negative(x)) fprintf(f, "  -");
        else  fprintf(f, "  0");
    }
    fprintf(f,"\n");
  }
  fprintf(f,"end\n");
  dd_clear(x);
}

void dd_WriteSignTableau2(FILE *f,dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
  dd_colindex nbindex_ref, dd_colindex nbindex,dd_rowindex bflag)
/* Write the sign tableau  A.T  and the reference basis */
{
  dd_colrange j;
  dd_rowrange i;
  mytype x;
  
  dd_init(x);
  fprintf(f," %ld  %ld  real\n",m_size,d_size);
  fprintf(f,"          |");
  for (j=1; j<= d_size; j++) fprintf(f,"%3ld",nbindex_ref[j]);
  fprintf(f,"\n          |");
  for (j=1; j<= d_size; j++) {
    fprintf(f,"%3ld",nbindex[j]);
  } fprintf(f,"\n  ------- | ");
  for (j=1; j<= d_size; j++) {
    fprintf(f,"---");
  } fprintf(f,"\n");
  for (i=1; i<= m_size; i++) {
    fprintf(f," %3ld(%3ld) |",i,bflag[i]);
    for (j=1; j<= d_size; j++) {
      dd_TableauEntry(&x,m_size,d_size,A,T,i,j);
      if (dd_Positive(x)) fprintf(f, "  +");
      else if (dd_Negative(x)) fprintf(f, "  -");
        else  fprintf(f, "  0");
    }
    fprintf(f,"\n");
  }
  fprintf(f,"end\n");
  dd_clear(x);
}


void dd_GetRedundancyInformation(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
  dd_colindex nbindex,dd_rowindex bflag, dd_rowset redset)
/* Some basic variables that are forced to be nonnegative will be output.  These are
   variables whose dictionary row components are all nonnegative.   */
{
  dd_colrange j;
  dd_rowrange i;
  mytype x;
  dd_boolean red=dd_FALSE,localdebug=dd_FALSE;
  long numbred=0;
  
  dd_init(x);
  for (i=1; i<= m_size; i++) {
    red=dd_TRUE;
    for (j=1; j<= d_size; j++) {
      dd_TableauEntry(&x,m_size,d_size,A,T,i,j);
      if (red && dd_Negative(x)) red=dd_FALSE;
    }
    if (bflag[i]<0 && red) {
      numbred+=1;
      set_addelem(redset,i);
    }
  }
  if (localdebug) fprintf(stderr,"\ndd_GetRedundancyInformation: %ld redundant rows over %ld\n",numbred, m_size);  
  dd_clear(x);
}


void dd_SelectDualSimplexPivot(dd_rowrange m_size,dd_colrange d_size,
    int Phase1,dd_Amatrix A,dd_Bmatrix T,dd_rowindex OV,
    dd_colindex nbindex_ref, dd_colindex nbindex,dd_rowindex bflag,
    dd_rowrange objrow,dd_colrange rhscol, dd_boolean lexicopivot,
    dd_rowrange *r,dd_colrange *s,int *selected,dd_LPStatusType *lps)
{ 
  /* selects a dual simplex pivot (*r,*s) if the current
     basis is dual feasible and not optimal. If not dual feasible,
     the procedure returns *selected=dd_FALSE and *lps=LPSundecided.
     If Phase1=dd_TRUE, the RHS column will be considered as the negative
     of the column of the largest variable (==m_size).  For this case, it is assumed
     that the caller used the auxiliary row (with variable m_size) to make the current
     dictionary dual feasible before calling this routine so that the nonbasic
     column for m_size corresponds to the auxiliary variable.
  */
  dd_boolean colselected=dd_FALSE,rowselected=dd_FALSE,
    dualfeasible=dd_TRUE,localdebug=dd_FALSE;
  dd_rowrange i,iref;
  dd_colrange j,k;
  mytype val,valn, minval,rat,minrat;
  static dd_Arow rcost;
  static dd_colrange d_last=0;
  static dd_colset tieset,stieset;  /* store the column indices with tie */

  dd_init(val); dd_init(valn); dd_init(minval); dd_init(rat); dd_init(minrat);
  if (d_last<d_size) {
    if (d_last>0) {
      for (j=1; j<=d_last; j++){ dd_clear(rcost[j-1]);}
      free(rcost);
      set_free(tieset);
      set_free(stieset);
    }
    rcost=(mytype*) calloc(d_size,sizeof(mytype));
    for (j=1; j<=d_size; j++){ dd_init(rcost[j-1]);}
    set_initialize(&tieset,d_size);
    set_initialize(&stieset,d_size);
  }
  d_last=d_size;

  *r=0; *s=0;
  *selected=dd_FALSE;
  *lps=dd_LPSundecided;
  for (j=1; j<=d_size; j++){
    if (j!=rhscol){
      dd_TableauEntry(&(rcost[j-1]),m_size,d_size,A,T,objrow,j);
      if (dd_Positive(rcost[j-1])) { 
        dualfeasible=dd_FALSE;
      }
    }
  }
  if (dualfeasible){
    while ((*lps==dd_LPSundecided) && (!rowselected) && (!colselected)) {
      for (i=1; i<=m_size; i++) {
        if (i!=objrow && bflag[i]==-1) {  /* i is a basic variable */
          if (Phase1){
            dd_TableauEntry(&val, m_size,d_size,A,T,i,bflag[m_size]);
            dd_neg(val,val);
            /* for dual Phase I.  The RHS (dual objective) is the negative of the auxiliary variable column. */
          } 
          else {dd_TableauEntry(&val,m_size,d_size,A,T,i,rhscol);}
          if (dd_Smaller(val,minval)) {
            *r=i;
            dd_set(minval,val);
          }
        }
      }
      if (dd_Nonnegative(minval)) {
        *lps=dd_Optimal;
      }
      else {
        rowselected=dd_TRUE;
        set_emptyset(tieset);
        for (j=1; j<=d_size; j++){
          dd_TableauEntry(&val,m_size,d_size,A,T,*r,j);
          if (j!=rhscol && dd_Positive(val)) {
            dd_div(rat,rcost[j-1],val);
            dd_neg(rat,rat);
            if (*s==0 || dd_Smaller(rat,minrat)){
              dd_set(minrat,rat);
              *s=j;
              set_emptyset(tieset);
              set_addelem(tieset, j);
            } else if (dd_Equal(rat,minrat)){
              set_addelem(tieset,j);
            }
          }
        }
        if (*s>0) {
          if (!lexicopivot || set_card(tieset)==1){
            colselected=dd_TRUE; *selected=dd_TRUE;
          } else { /* lexicographic rule with respect to the given reference cobasis.  */
            if (localdebug) {printf("Tie occurred at:"); set_write(tieset); printf("\n");
              dd_WriteTableau(stderr,m_size,d_size,A,T,nbindex,bflag);
            }
            *s=0;
            k=2; /* k runs through the column indices except RHS. */
            do {
              iref=nbindex_ref[k];  /* iref runs though the reference basic indices */
              if (iref>0) {
                j=bflag[iref];
                if (j>0) {
                  if (set_member(j,tieset) && set_card(tieset)==1) {
                    *s=j;
                     colselected=dd_TRUE;
                  } else {
                    set_delelem(tieset, j);
                    /* iref is cobasic, and the corresponding col is not the pivot column except it is the last one. */
                  }
                } else {
                  *s=0;
                  for (j=1; j<=d_size; j++){
                    if (set_member(j,tieset)) {
                      dd_TableauEntry(&val,m_size,d_size,A,T,*r,j);
                      dd_TableauEntry(&valn,m_size,d_size,A,T,iref,j);
                      if (j!=rhscol && dd_Positive(val)) {
                        dd_div(rat,valn,val);
                        if (*s==0 || dd_Smaller(rat,minrat)){
                          dd_set(minrat,rat);
                          *s=j;
                          set_emptyset(stieset);
                          set_addelem(stieset, j);
                        } else if (dd_Equal(rat,minrat)){
                          set_addelem(stieset,j);
                        }
                      }
                    }
                  }
                  set_copy(tieset,stieset);              
                  if (set_card(tieset)==1) colselected=dd_TRUE;
                }
              }
              k+=1;
            } while (!colselected && k<=d_size);
            *selected=dd_TRUE;
          }
        } else *lps=dd_Inconsistent;
      }
    } /* end of while */
  }
  if (localdebug) {
     if (Phase1) fprintf(stderr,"Phase 1 : select %ld,%ld\n",*r,*s);
     else fprintf(stderr,"Phase 2 : select %ld,%ld\n",*r,*s);
  }
  dd_clear(val); dd_clear(valn); dd_clear(minval); dd_clear(rat); dd_clear(minrat);
}

void dd_TableauEntry(mytype *x,dd_rowrange m_size, dd_colrange d_size, dd_Amatrix X, dd_Bmatrix T,
				dd_rowrange r, dd_colrange s)
/* Compute the (r,s) entry of X.T   */
{
  dd_colrange j;
  mytype temp;

  dd_init(temp);
  dd_set(*x,dd_purezero);
  for (j=0; j< d_size; j++) {
    dd_mul(temp,X[r-1][j], T[j][s-1]);
    dd_add(*x, *x, temp);
  }
  dd_clear(temp);
}

void dd_SelectPivot2(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,dd_Bmatrix T,
            dd_RowOrderType roworder,dd_rowindex ordervec, rowset equalityset,
            dd_rowrange rowmax,rowset NopivotRow,
            colset NopivotCol,dd_rowrange *r,dd_colrange *s,
            dd_boolean *selected)
/* Select a position (*r,*s) in the matrix A.T such that (A.T)[*r][*s] is nonzero
   The choice is feasible, i.e., not on NopivotRow and NopivotCol, and
   best with respect to the specified roworder 
 */
{
  int stop;
  dd_rowrange i,rtemp;
  rowset rowexcluded;
  mytype Xtemp;
  dd_boolean localdebug=dd_FALSE;

  stop = dd_FALSE;
  localdebug=dd_debug;
  dd_init(Xtemp);
  set_initialize(&rowexcluded,m_size);
  set_copy(rowexcluded,NopivotRow);
  for (i=rowmax+1;i<=m_size;i++) {
    set_addelem(rowexcluded,i);   /* cannot pivot on any row > rmax */
  }
  *selected = dd_FALSE;
  do {
    rtemp=0; i=1;
    while (i<=m_size && rtemp==0) {  /* equalityset vars have highest priorities */
      if (set_member(i,equalityset) && !set_member(i,rowexcluded)){
        if (localdebug) fprintf(stderr,"marked set %ld chosen as a candidate\n",i);
        rtemp=i;
      }
      i++;
    }
    if (rtemp==0) dd_SelectPreorderedNext2(m_size,d_size,rowexcluded,ordervec,&rtemp);;
    if (rtemp>=1) {
      *r=rtemp;
      *s=1;
      while (*s <= d_size && !*selected) {
        dd_TableauEntry(&Xtemp,m_size,d_size,A,T,*r,*s);
        if (!set_member(*s,NopivotCol) && dd_Nonzero(Xtemp)) {
          *selected = dd_TRUE;
          stop = dd_TRUE;
        } else {
          (*s)++;
        }
      }
      if (!*selected) {
        set_addelem(rowexcluded,rtemp);
      }
    }
    else {
      *r = 0;
      *s = 0;
      stop = dd_TRUE;
    }
  } while (!stop);
  set_free(rowexcluded); dd_clear(Xtemp);
}

void dd_GaussianColumnPivot(dd_rowrange m_size, dd_colrange d_size, 
    dd_Amatrix X, dd_Bmatrix T, dd_rowrange r, dd_colrange s)
/* Update the Transformation matrix T with the pivot operation on (r,s) 
   This procedure performs a implicit pivot operation on the matrix X by
   updating the dual basis inverse  T.
 */
{
  dd_colrange j, j1;
  mytype Xtemp0, Xtemp1, Xtemp;
  static dd_Arow Rtemp;
  static dd_colrange last_d=0;

  dd_init(Xtemp0); dd_init(Xtemp1); dd_init(Xtemp);
  if (last_d!=d_size){
    if (last_d>0) {
      for (j=1; j<=last_d; j++) dd_clear(Rtemp[j-1]);
      free(Rtemp);
    }
    Rtemp=(mytype*)calloc(d_size,sizeof(mytype));
    for (j=1; j<=d_size; j++) dd_init(Rtemp[j-1]);
    last_d=d_size;
  }

  for (j=1; j<=d_size; j++) {
    dd_TableauEntry(&(Rtemp[j-1]), m_size, d_size, X, T, r,j);
  }
  dd_set(Xtemp0,Rtemp[s-1]);
  for (j = 1; j <= d_size; j++) {
    if (j != s) {
      dd_div(Xtemp,Rtemp[j-1],Xtemp0);
      dd_set(Xtemp1,dd_purezero);
      for (j1 = 1; j1 <= d_size; j1++){
        dd_mul(Xtemp1,Xtemp,T[j1-1][s - 1]);
        dd_sub(T[j1-1][j-1],T[j1-1][j-1],Xtemp1);
 /*     T[j1-1][j-1] -= T[j1-1][s - 1] * Xtemp / Xtemp0;  */
      }
    }
  }
  for (j = 1; j <= d_size; j++)
    dd_div(T[j-1][s - 1],T[j-1][s - 1],Xtemp0);

  dd_clear(Xtemp0); dd_clear(Xtemp1); dd_clear(Xtemp);
}

void dd_GaussianColumnPivot2(dd_rowrange m_size,dd_colrange d_size,
    dd_Amatrix A,dd_Bmatrix T,dd_colindex nbindex,dd_rowindex bflag,dd_rowrange r,dd_colrange s)
/* Update the Transformation matrix T with the pivot operation on (r,s) 
   This procedure performs a implicit pivot operation on the matrix A by
   updating the dual basis inverse  T.
 */
{
  int localdebug=dd_FALSE;
  long entering;

  if (dd_debug) localdebug=dd_debug;