Added cddlib-094h from http://www.inf.ethz.ch/personal/fukudak/cdd_home/
Change-Id: I64519509269e434b1b9ea87c3fe0805e711c0ac9
diff --git a/third_party/cddlib/lib-src/cddproj.c b/third_party/cddlib/lib-src/cddproj.c
new file mode 100644
index 0000000..2ef7e7a
--- /dev/null
+++ b/third_party/cddlib/lib-src/cddproj.c
@@ -0,0 +1,235 @@
+/* cddproj.c: Polyhedral Projections in cddlib
+ written by Komei Fukuda, fukuda@cs.mcgill.ca
+ Version 0.94, Aug. 4, 2005
+*/
+
+/* cddlib : C-library of the double description method for
+ computing all vertices and extreme rays of the polyhedron
+ P= {x : b - A x >= 0}.
+ Please read COPYING (GNU General Public Licence) and
+ the manual cddlibman.tex for detail.
+*/
+
+#include "setoper.h" /* set operation library header (Ver. June 1, 2000 or later) */
+#include "cdd.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <math.h>
+#include <string.h>
+
+dd_MatrixPtr dd_BlockElimination(dd_MatrixPtr M, dd_colset delset, dd_ErrorType *error)
+/* Eliminate the variables (columns) delset by
+ the Block Elimination with dd_DoubleDescription algorithm.
+
+ Given (where y is to be eliminated):
+ c1 + A1 x + B1 y >= 0
+ c2 + A2 x + B2 y = 0
+
+ 1. First construct the dual system: z1^T B1 + z2^T B2 = 0, z1 >= 0.
+ 2. Compute the generators of the dual.
+ 3. Then take the linear combination of the original system with each generator.
+ 4. Remove redundant inequalies.
+
+*/
+{
+ dd_MatrixPtr Mdual=NULL, Mproj=NULL, Gdual=NULL;
+ dd_rowrange i,h,m,mproj,mdual,linsize;
+ dd_colrange j,k,d,dproj,ddual,delsize;
+ dd_colindex delindex;
+ mytype temp,prod;
+ dd_PolyhedraPtr dualpoly;
+ dd_ErrorType err=dd_NoError;
+ dd_boolean localdebug=dd_FALSE;
+
+ *error=dd_NoError;
+ m= M->rowsize;
+ d= M->colsize;
+ delindex=(long*)calloc(d+1,sizeof(long));
+ dd_init(temp);
+ dd_init(prod);
+
+ k=0; delsize=0;
+ for (j=1; j<=d; j++){
+ if (set_member(j, delset)){
+ k++; delsize++;
+ delindex[k]=j; /* stores the kth deletion column index */
+ }
+ }
+ if (localdebug) dd_WriteMatrix(stdout, M);
+
+ linsize=set_card(M->linset);
+ ddual=m+1;
+ mdual=delsize + m - linsize; /* #equalitions + dimension of z1 */
+
+ /* setup the dual matrix */
+ Mdual=dd_CreateMatrix(mdual, ddual);
+ Mdual->representation=dd_Inequality;
+ for (i = 1; i <= delsize; i++){
+ set_addelem(Mdual->linset,i); /* equality */
+ for (j = 1; j <= m; j++) {
+ dd_set(Mdual->matrix[i-1][j], M->matrix[j-1][delindex[i]-1]);
+ }
+ }
+
+ k=0;
+ for (i = 1; i <= m; i++){
+ if (!set_member(i, M->linset)){
+ /* set nonnegativity for the dual variable associated with
+ each non-linearity inequality. */
+ k++;
+ dd_set(Mdual->matrix[delsize+k-1][i], dd_one);
+ }
+ }
+
+ /* 2. Compute the generators of the dual system. */
+ dualpoly=dd_DDMatrix2Poly(Mdual, &err);
+ Gdual=dd_CopyGenerators(dualpoly);
+
+ /* 3. Take the linear combination of the original system with each generator. */
+ dproj=d-delsize;
+ mproj=Gdual->rowsize;
+ Mproj=dd_CreateMatrix(mproj, dproj);
+ Mproj->representation=dd_Inequality;
+ set_copy(Mproj->linset, Gdual->linset);
+
+ for (i=1; i<=mproj; i++){
+ k=0;
+ for (j=1; j<=d; j++){
+ if (!set_member(j, delset)){
+ k++; /* new index of the variable x_j */
+ dd_set(prod, dd_purezero);
+ for (h = 1; h <= m; h++){
+ dd_mul(temp,M->matrix[h-1][j-1],Gdual->matrix[i-1][h]);
+ dd_add(prod,prod,temp);
+ }
+ dd_set(Mproj->matrix[i-1][k-1],prod);
+ }
+ }
+ }
+ if (localdebug) printf("Size of the projection system: %ld x %ld\n", mproj, dproj);
+
+ dd_FreePolyhedra(dualpoly);
+ free(delindex);
+ dd_clear(temp);
+ dd_clear(prod);
+ dd_FreeMatrix(Mdual);
+ dd_FreeMatrix(Gdual);
+ return Mproj;
+}
+
+
+dd_MatrixPtr dd_FourierElimination(dd_MatrixPtr M,dd_ErrorType *error)
+/* Eliminate the last variable (column) from the given H-matrix using
+ the standard Fourier Elimination.
+ */
+{
+ dd_MatrixPtr Mnew=NULL;
+ dd_rowrange i,inew,ip,in,iz,m,mpos=0,mneg=0,mzero=0,mnew;
+ dd_colrange j,d,dnew;
+ dd_rowindex posrowindex, negrowindex,zerorowindex;
+ mytype temp1,temp2;
+ dd_boolean localdebug=dd_FALSE;
+
+ *error=dd_NoError;
+ m= M->rowsize;
+ d= M->colsize;
+ if (d<=1){
+ *error=dd_ColIndexOutOfRange;
+ if (localdebug) {
+ printf("The number of column is too small: %ld for Fourier's Elimination.\n",d);
+ }
+ goto _L99;
+ }
+
+ if (M->representation==dd_Generator){
+ *error=dd_NotAvailForV;
+ if (localdebug) {
+ printf("Fourier's Elimination cannot be applied to a V-polyhedron.\n");
+ }
+ goto _L99;
+ }
+
+ if (set_card(M->linset)>0){
+ *error=dd_CannotHandleLinearity;
+ if (localdebug) {
+ printf("The Fourier Elimination function does not handle equality in this version.\n");
+ }
+ goto _L99;
+ }
+
+ /* Create temporary spaces to be removed at the end of this function */
+ posrowindex=(long*)calloc(m+1,sizeof(long));
+ negrowindex=(long*)calloc(m+1,sizeof(long));
+ zerorowindex=(long*)calloc(m+1,sizeof(long));
+ dd_init(temp1);
+ dd_init(temp2);
+
+ for (i = 1; i <= m; i++) {
+ if (dd_Positive(M->matrix[i-1][d-1])){
+ mpos++;
+ posrowindex[mpos]=i;
+ } else if (dd_Negative(M->matrix[i-1][d-1])) {
+ mneg++;
+ negrowindex[mneg]=i;
+ } else {
+ mzero++;
+ zerorowindex[mzero]=i;
+ }
+ } /*of i*/
+
+ if (localdebug) {
+ dd_WriteMatrix(stdout, M);
+ printf("No of (+ - 0) rows = (%ld, %ld, %ld)\n", mpos,mneg, mzero);
+ }
+
+ /* The present code generates so many redundant inequalities and thus
+ is quite useless, except for very small examples
+ */
+ mnew=mzero+mpos*mneg; /* the total number of rows after elimination */
+ dnew=d-1;
+
+ Mnew=dd_CreateMatrix(mnew, dnew);
+ dd_CopyArow(Mnew->rowvec, M->rowvec, dnew);
+/* set_copy(Mnew->linset,M->linset); */
+ Mnew->numbtype=M->numbtype;
+ Mnew->representation=M->representation;
+ Mnew->objective=M->objective;
+
+
+ /* Copy the inequalities independent of x_d to the top of the new matrix. */
+ for (iz = 1; iz <= mzero; iz++){
+ for (j = 1; j <= dnew; j++) {
+ dd_set(Mnew->matrix[iz-1][j-1], M->matrix[zerorowindex[iz]-1][j-1]);
+ }
+ }
+
+ /* Create the new inequalities by combining x_d positive and negative ones. */
+ inew=mzero; /* the index of the last x_d zero inequality */
+ for (ip = 1; ip <= mpos; ip++){
+ for (in = 1; in <= mneg; in++){
+ inew++;
+ dd_neg(temp1, M->matrix[negrowindex[in]-1][d-1]);
+ for (j = 1; j <= dnew; j++) {
+ dd_LinearComb(temp2,M->matrix[posrowindex[ip]-1][j-1],temp1,\
+ M->matrix[negrowindex[in]-1][j-1],\
+ M->matrix[posrowindex[ip]-1][d-1]);
+ dd_set(Mnew->matrix[inew-1][j-1],temp2);
+ }
+ dd_Normalize(dnew,Mnew->matrix[inew-1]);
+ }
+ }
+
+
+ free(posrowindex);
+ free(negrowindex);
+ free(zerorowindex);
+ dd_clear(temp1);
+ dd_clear(temp2);
+
+ _L99:
+ return Mnew;
+}
+
+
+/* end of cddproj.c */