Sandia Home Sandia Home
Main Page | Publications | Downloads | Configuration | Running the Code | Solver Parameters | FAQ | Namespace List | Class Hierarchy | Class List | File List | Namespace Members | Class Members | File Members | Related Pages

cddproj.c

Go to the documentation of this file.
00001 
00043 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00044 
00045 /* cddproj.c:  Polyhedral Projections in cddlib
00046    written by Komei Fukuda, fukuda@cs.mcgill.ca
00047    Version 0.94, Aug. 4, 2005
00048 */
00049 
00050 /* cddlib : C-library of the double description method for
00051    computing all vertices and extreme rays of the polyhedron 
00052    P= {x :  b - A x >= 0}.  
00053    Please read COPYING (GNU General Public Licence) and
00054    the manual cddlibman.tex for detail.
00055 */
00056 
00057 #include "setoper.h"  /* set operation library header (Ver. June 1, 2000 or later) */
00058 #include "cdd.h"
00059 #include <stdio.h>
00060 #include <stdlib.h>
00061 #include <time.h>
00062 #include <math.h>
00063 #include <string.h>
00064 
00065 dd_MatrixPtr dd_BlockElimination(dd_MatrixPtr M, dd_colset delset, dd_ErrorType *error)
00066 /* Eliminate the variables (columns) delset by
00067    the Block Elimination with dd_DoubleDescription algorithm.
00068 
00069    Given (where y is to be eliminated):
00070    c1 + A1 x + B1 y >= 0
00071    c2 + A2 x + B2 y =  0
00072 
00073    1. First construct the dual system:  z1^T B1 + z2^T B2 = 0, z1 >= 0.
00074    2. Compute the generators of the dual.
00075    3. Then take the linear combination of the original system with each generator.
00076    4. Remove redundant inequalies.
00077 
00078 */
00079 {
00080   dd_MatrixPtr Mdual=NULL, Mproj=NULL, Gdual=NULL;
00081   dd_rowrange i,h,m,mproj,mdual,linsize;
00082   dd_colrange j,k,d,dproj,ddual,delsize;
00083   dd_colindex delindex;
00084   mytype temp,prod;
00085   dd_PolyhedraPtr dualpoly;
00086   dd_ErrorType err=dd_NoError;
00087   dd_boolean localdebug=dd_FALSE;
00088 
00089   *error=dd_NoError;
00090   m= M->rowsize;
00091   d= M->colsize;
00092   delindex=(long*)calloc(d+1,sizeof(long));
00093   dd_init(temp);
00094   dd_init(prod);
00095 
00096   k=0; delsize=0;
00097   for (j=1; j<=d; j++){
00098     if (set_member(j, delset)){
00099       k++;  delsize++;
00100       delindex[k]=j;  /* stores the kth deletion column index */
00101     }
00102   }
00103   if (localdebug) dd_WriteMatrix(stdout, M);
00104 
00105   linsize=set_card(M->linset);
00106   ddual=m+1;
00107   mdual=delsize + m - linsize;  /* #equalitions + dimension of z1 */
00108 
00109   /* setup the dual matrix */
00110   Mdual=dd_CreateMatrix(mdual, ddual);
00111   Mdual->representation=dd_Inequality;
00112   for (i = 1; i <= delsize; i++){
00113     set_addelem(Mdual->linset,i);  /* equality */
00114     for (j = 1; j <= m; j++) {
00115       dd_set(Mdual->matrix[i-1][j], M->matrix[j-1][delindex[i]-1]);
00116     }
00117   } 
00118 
00119   k=0;
00120   for (i = 1; i <= m; i++){
00121     if (!set_member(i, M->linset)){
00122       /* set nonnegativity for the dual variable associated with
00123          each non-linearity inequality. */
00124       k++;
00125       dd_set(Mdual->matrix[delsize+k-1][i], dd_one);  
00126     }
00127   } 
00128   
00129   /* 2. Compute the generators of the dual system. */
00130   dualpoly=dd_DDMatrix2Poly(Mdual, &err);
00131   Gdual=dd_CopyGenerators(dualpoly);
00132 
00133   /* 3. Take the linear combination of the original system with each generator.  */
00134   dproj=d-delsize;
00135   mproj=Gdual->rowsize;
00136   Mproj=dd_CreateMatrix(mproj, dproj);
00137   Mproj->representation=dd_Inequality;
00138   set_copy(Mproj->linset, Gdual->linset);
00139 
00140   for (i=1; i<=mproj; i++){
00141     k=0;
00142     for (j=1; j<=d; j++){
00143       if (!set_member(j, delset)){
00144         k++;  /* new index of the variable x_j  */
00145         dd_set(prod, dd_purezero);
00146         for (h = 1; h <= m; h++){
00147           dd_mul(temp,M->matrix[h-1][j-1],Gdual->matrix[i-1][h]); 
00148           dd_add(prod,prod,temp);
00149         }
00150         dd_set(Mproj->matrix[i-1][k-1],prod);
00151       }
00152     }
00153   }
00154   if (localdebug) printf("Size of the projection system: %ld x %ld\n", mproj, dproj);
00155   
00156   dd_FreePolyhedra(dualpoly);
00157   free(delindex);
00158   dd_clear(temp);
00159   dd_clear(prod);
00160   dd_FreeMatrix(Mdual);
00161   dd_FreeMatrix(Gdual);
00162   return Mproj;
00163 }
00164 
00165 
00166 dd_MatrixPtr dd_FourierElimination(dd_MatrixPtr M,dd_ErrorType *error)
00167 /* Eliminate the last variable (column) from the given H-matrix using 
00168    the standard Fourier Elimination.
00169  */
00170 {
00171   dd_MatrixPtr Mnew=NULL;
00172   dd_rowrange i,inew,ip,in,iz,m,mpos=0,mneg=0,mzero=0,mnew;
00173   dd_colrange j,d,dnew;
00174   dd_rowindex posrowindex, negrowindex,zerorowindex;
00175   mytype temp1,temp2;
00176   dd_boolean localdebug=dd_FALSE;
00177 
00178   *error=dd_NoError;
00179   m= M->rowsize;
00180   d= M->colsize;
00181   if (d<=1){
00182     *error=dd_ColIndexOutOfRange;
00183     if (localdebug) {
00184       printf("The number of column is too small: %ld for Fourier's Elimination.\n",d);
00185     }
00186     goto _L99;
00187   }
00188 
00189   if (M->representation==dd_Generator){
00190     *error=dd_NotAvailForV;
00191     if (localdebug) {
00192       printf("Fourier's Elimination cannot be applied to a V-polyhedron.\n");
00193     }
00194     goto _L99;
00195   }
00196 
00197   if (set_card(M->linset)>0){
00198     *error=dd_CannotHandleLinearity;
00199     if (localdebug) {
00200       printf("The Fourier Elimination function does not handle equality in this version.\n");
00201     }
00202     goto _L99;
00203   }
00204 
00205   /* Create temporary spaces to be removed at the end of this function */
00206   posrowindex=(long*)calloc(m+1,sizeof(long));
00207   negrowindex=(long*)calloc(m+1,sizeof(long));
00208   zerorowindex=(long*)calloc(m+1,sizeof(long));
00209   dd_init(temp1);
00210   dd_init(temp2);
00211 
00212   for (i = 1; i <= m; i++) {
00213     if (dd_Positive(M->matrix[i-1][d-1])){
00214       mpos++;
00215       posrowindex[mpos]=i;
00216     } else if (dd_Negative(M->matrix[i-1][d-1])) {
00217       mneg++;
00218       negrowindex[mneg]=i;
00219     } else {
00220       mzero++;
00221       zerorowindex[mzero]=i;
00222     }
00223   }  /*of i*/
00224 
00225   if (localdebug) {
00226     dd_WriteMatrix(stdout, M);
00227     printf("No of  (+  -  0) rows = (%ld, %ld, %ld)\n", mpos,mneg, mzero);
00228   }
00229 
00230   /* The present code generates so many redundant inequalities and thus
00231      is quite useless, except for very small examples
00232   */
00233   mnew=mzero+mpos*mneg;  /* the total number of rows after elimination */
00234   dnew=d-1;
00235 
00236   Mnew=dd_CreateMatrix(mnew, dnew);
00237   dd_CopyArow(Mnew->rowvec, M->rowvec, dnew);
00238 /*  set_copy(Mnew->linset,M->linset);  */
00239   Mnew->numbtype=M->numbtype;
00240   Mnew->representation=M->representation;
00241   Mnew->objective=M->objective;
00242 
00243 
00244   /* Copy the inequalities independent of x_d to the top of the new matrix. */
00245   for (iz = 1; iz <= mzero; iz++){
00246     for (j = 1; j <= dnew; j++) {
00247       dd_set(Mnew->matrix[iz-1][j-1], M->matrix[zerorowindex[iz]-1][j-1]);
00248     }
00249   } 
00250 
00251   /* Create the new inequalities by combining x_d positive and negative ones. */
00252   inew=mzero;  /* the index of the last x_d zero inequality */
00253   for (ip = 1; ip <= mpos; ip++){
00254     for (in = 1; in <= mneg; in++){
00255       inew++;
00256       dd_neg(temp1, M->matrix[negrowindex[in]-1][d-1]);
00257       for (j = 1; j <= dnew; j++) {
00258         dd_LinearComb(temp2,M->matrix[posrowindex[ip]-1][j-1],temp1,\
00259           M->matrix[negrowindex[in]-1][j-1],\
00260           M->matrix[posrowindex[ip]-1][d-1]);
00261         dd_set(Mnew->matrix[inew-1][j-1],temp2);
00262       }
00263       dd_Normalize(dnew,Mnew->matrix[inew-1]);
00264     }
00265   } 
00266 
00267 
00268   free(posrowindex);
00269   free(negrowindex);
00270   free(zerorowindex);
00271   dd_clear(temp1);
00272   dd_clear(temp2);
00273 
00274  _L99:
00275   return Mnew;
00276 }
00277 
00278 
00279 /* end of cddproj.c */
00280 
00281 #endif

 

© Sandia Corporation | Site Contact | Privacy and Security

Generated on Fri Feb 16 10:33:35 2007 for APPSPACK 5.0.1 by doxygen 1.3.9.1 written by Dimitri van Heesch, © 1997-2002