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

APPSPACK_CDDLIB_Interface.c

Go to the documentation of this file.
00001 
00002 /*
00003 //@HEADER
00004 // ************************************************************************
00005 // 
00006 //          APPSPACK: Asynchronous Parallel Pattern Search
00007 //                 Copyright (2003) Sandia Corporation
00008 // 
00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00010 // license for use of this work by or on behalf of the U.S. Government.
00011 // 
00012 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //                                                                                 
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA.                                                                           .
00026 // 
00027 // Questions? Contact Tammy Kolda (tgkolda@sandia.gov) 
00028 // 
00029 // ************************************************************************
00030 //@HEADER
00031 */
00032 
00033 #include "APPSPACK_CDDLIB_Interface.h"
00034 
00035 #ifdef HAVE_CDDLIB
00036 
00043 /* Includes from cddlib. */
00044 #include "setoper.h"
00045 #include "cdd.h"
00046 
00047 #ifndef  DOXYGEN_SHOULD_SKIP_THIS
00048 #define OK 0
00049 #define OKEY_DOKEY  0
00050 #define ACK        16
00051 #endif
00052 
00053 int compute_cone_generators(int *num_pointy,    double ***P,
00054                             int *num_lineality, double ***L,
00055                             int n, 
00056                             int num_equalities,   double **Eq,
00057                             int num_inequalities, double **Iq,
00058                             int append)
00059 {
00060   static dd_PolyhedraPtr poly = NULL;
00061   static int last_n;
00062 
00063   dd_MatrixPtr A, G;
00064   dd_ErrorType err;
00065 
00066   int i, j, k;
00067   int linsize;
00068   int nl, np;
00069 
00070   int num_constraints;
00071 
00072   /* Some local initializations. */
00073 
00074   *P = NULL;  *num_pointy = 0;
00075   *L = NULL;  *num_lineality = 0;
00076 
00077   num_constraints = num_equalities + num_inequalities;
00078 
00079   /* Sanity check. */
00080 
00081   if (num_constraints <= 0) {
00082     printf("Bogus number of constraints specified: %d !!\n", num_constraints);
00083     return (ACK);
00084   }
00085 
00086   /* Initialize cddlib. */
00087 
00088   if (poly == NULL) {
00089     dd_set_global_constants();
00090   }
00091 
00092   if (append) {
00093     if (poly == NULL) {
00094       printf("%s: 'append' makes no sense---this is the first call\n",
00095              __FILE__);
00096       return (ACK);
00097     }
00098     if (n != last_n) {
00099       printf("%s: 'append' makes no sense---the dimension has changed\n",
00100              __FILE__);
00101       return (ACK);
00102     }
00103   }
00104 
00105   /* 
00106    * If we've called cddlib before and we are starting from scratch with a 
00107    * new set of constraints, then free the state we have retained from the
00108    * previous call.
00109    */
00110 
00111   if (!append && (poly != NULL)) {
00112     dd_FreePolyhedra(poly);
00113   }
00114 
00115   /*
00116    * Slog through the computation of the generators of T(x,r)...
00117    *
00118    * cddlib expects the constraints to have the form 
00119    *    Ax <= b   <=>   b - Ax >= 0,
00120    * and expects to be given b and -A.  Note the sign reversal of the 
00121    * constraint coefficients!  Tricky, tricky, tricky!
00122    *
00123    * The RHS b goes in A->matrix(:,0) while -A goes in A->matrix(:,1:n+1).
00124    */
00125 
00126   A = dd_CreateMatrix(num_constraints, n+1);
00127   A->representation = dd_Inequality;
00128 
00129   for (i = 0; i < num_equalities; i++) {
00130     dd_set_si(A->matrix[i][0], 0);  /* For our problem, b = 0. */
00131     for (j = 0; j < n; j++) {
00132       dd_set_d(A->matrix[i][j+1], -Eq[i][j]);
00133     }
00134     set_addelem(A->linset, i+1);
00135   }
00136   for (k = 0; i < num_constraints; i++, k++) {
00137     dd_set_si(A->matrix[i][0], 0);  /* For our problem, b = 0. */
00138     for (j = 0; j < n; j++) {
00139       dd_set_d(A->matrix[i][j+1], -Iq[k][j]);
00140     }
00141   }
00142 
00143   if (!append) {
00144     poly = dd_DDMatrix2Poly(A, &err);
00145   }
00146   else {
00147     dd_DDInputAppend(&poly, A, &err);
00148   }
00149 
00150   if (err != dd_NoError) {
00151     return (ACK);
00152   }
00153 
00154   /* Translate the generator representation to a matrix format. */
00155   G = dd_CopyGenerators(poly);
00156 
00157   if (G == NULL) {
00158     printf("%s: G == NULL encountered at line %d\n", __FILE__, __LINE__);
00159     return (ACK);
00160   }
00161 
00162   if ((G->colsize) != (n+1)) { /* This should never fail... */
00163     printf("%s: G->colsize != n+1 encountered at line %d\n", 
00164            __FILE__, __LINE__);
00165     return (ACK);
00166   }
00167 
00168   /*
00169    * Allocate space for the generators.
00170    */
00171 
00172   /* linsize is the dimension of the lineality space. */
00173 
00174   linsize = set_card(G->linset);
00175 
00176   *num_pointy    = G->rowsize - linsize;
00177   *num_lineality = linsize;
00178 
00179   /* Watch out!  If the input generates all of R^{n}, then what */
00180   /* appears in the pointy part is just the vertex at 0!! */
00181   printf("dim. of lineality space:            %d\n", *num_lineality);
00182   printf("no. of extreme rays in pointy part: %d\n", *num_pointy);
00183 
00184   if (*num_pointy > 0) {
00185     *P = (double**) malloc( (*num_pointy) * sizeof(double*) );
00186     if (*P == NULL) {
00187       printf("%s: Memory allocation failed! (line %d)\n", __FILE__, __LINE__);
00188       printf("%s: Requested %d bytes\n", __FILE__, 
00189              (*num_pointy) * sizeof(double*));
00190       return (ACK);
00191     }
00192   }
00193 
00194   if (*num_lineality > 0) {
00195     *L = (double**) malloc( (*num_lineality) * sizeof(double*) );
00196     if (*L == NULL) {
00197       printf("%s: Memory allocation failed! (line %d)\n",
00198              __FILE__, __LINE__);
00199       printf("%s: Requested %d bytes\n", __FILE__,
00200              (*num_lineality) * sizeof(double*));
00201       return (ACK);
00202     }
00203   }
00204 
00205   for (i = 0, np = 0, nl = 0; i < G->rowsize; i++) {
00206     
00207     /* 
00208      * Check whether this generator is part of the lineality space.
00209      */
00210 
00211     if (set_member(i+1, G->linset)) {
00212       (*L)[nl] = (double*) malloc(n*sizeof(double));
00213       if ((*L)[nl] == NULL) {
00214         printf("%s: Memory allocation failed! (line %d)\n",
00215                __FILE__, __LINE__);
00216         printf("%s: Requested %d bytes\n", __FILE__, n * sizeof(double*));
00217         return (ACK);
00218       }
00219       for (j = 0; j < n; j++) {
00220         (*L)[nl][j] = dd_get_d(G->matrix[i][j+1]);
00221       }
00222       nl++;
00223     }
00224     else {
00225       (*P)[np] = (double*) malloc(n*sizeof(double));
00226       if ((*P)[np] == NULL) {
00227         printf("%s: Memory allocation failed! (line %d)\n",
00228                __FILE__, __LINE__);
00229         printf("%s: Requested %d bytes\n", __FILE__, n * sizeof(double*));
00230         return (ACK);
00231       }
00232       for (j = 0; j < n; j++) {
00233         (*P)[np][j] = dd_get_d(G->matrix[i][j+1]);
00234       }
00235       np++;
00236     }
00237   }
00238 
00239   /* Clean up! */
00240 
00241   dd_FreeMatrix(A);
00242   dd_FreeMatrix(G);
00243 
00244   last_n = n;
00245   return (OKEY_DOKEY);
00246 }
00247 
00248 
00249 /*****************************************************************************
00250  * A test driver, included at no extra charge!!
00251  */
00252 
00253 #if defined(TEST_MAIN)
00254 #include "math.h"
00255 
00256 int main (int argc, char **argv)
00257 {
00258   int i, j, k;
00259 
00260   int retcd;
00261   int num_pointy, num_lineality;
00262   double **P;
00263   double **L;
00264   double **Eq, **Iq;
00265   sp47opts opts = {true, true, false, 1.0e-08};
00266 
00267   int
00268     num_equalities,
00269     num_inequalities;
00270   int n;
00271 
00272   num_equalities   = 3;
00273   num_inequalities = 0;
00274   n = 7;
00275 
00276   Eq = (double**) malloc(num_equalities*sizeof(double*));
00277   for (i = 0; i < num_equalities; i++) {
00278     Eq[i] = (double*) malloc(n*sizeof(double)); 
00279     for (j = 0; j < n; j++) {
00280       Eq[i][j] = drand48();
00281     }
00282   }
00283   Iq = NULL;
00284 
00285   for (k = 0; k < 100; k++) {
00286     compute_generators_new (&num_pointy, &P, &num_lineality, &L, 
00287                         n, num_equalities, Eq, num_inequalities, Iq, &opts);
00288 #if defined(VERBOSE)
00289     printf("Pointy directions:\n");
00290     for (i = 0; i < num_pointy; i++) {
00291       for (j = 0; j < n; j++) {
00292         printf("% f ", P[i][j]);
00293       }
00294       printf("\n");
00295     }
00296 
00297     printf("Lineality space:\n");
00298     for (i = 0; i < num_lineality; i++) {
00299       for (j = 0; j < n; j++) {
00300         printf("% f ", L[i][j]);
00301       }
00302       printf("\n");
00303     }
00304 #endif
00305 
00306     for (i = 0; i < num_pointy; i++) {
00307       free(P[i]);
00308     }
00309     free(P);
00310 
00311     for (i = 0; i < num_lineality; i++) {
00312       free(L[i]);
00313     }
00314     if (L != NULL) {
00315       free(L);
00316     }
00317   }
00318 
00319   for (i = 0; i < num_equalities; i++) {
00320     free(Eq[i]);
00321   }
00322   free(Eq);
00323 
00324   /***************************************************************
00325    * Test on a pointy cone.
00326    */
00327 
00328   num_equalities   = 0;
00329   num_inequalities = 20;
00330   n = 30;
00331 
00332   Eq = NULL;
00333   Iq = (double**) malloc(num_inequalities*sizeof(double*));
00334   for (i = 0; i < num_inequalities; i++) {
00335     Iq[i] = (double*) malloc(n*sizeof(double)); 
00336     for (j = 0; j < n; j++) {
00337       Iq[i][j] = drand48();
00338     }
00339   }
00340 
00341   for (k = 0; k < 500; k++) {
00342     compute_generators_new (&num_pointy, &P, &num_lineality, &L, 
00343                         n, num_equalities, Eq, num_inequalities, Iq, &opts);
00344 
00345 #if defined(VERBOSE)
00346     printf("Pointy directions:\n");
00347     for (i = 0; i < num_pointy; i++) {
00348       for (j = 0; j < n; j++) {
00349         printf("% f ", P[i][j]);
00350       }
00351       printf("\n");
00352     }
00353 
00354     printf("Lineality space:\n");
00355     for (i = 0; i < num_lineality; i++) {
00356       for (j = 0; j < n; j++) {
00357         printf("% f ", L[i][j]);
00358       }
00359       printf("\n");
00360     }
00361 #endif
00362 
00363     for (i = 0; i < num_pointy; i++) {
00364       free(P[i]);
00365     }
00366     free(P);
00367 
00368     for (i = 0; i < num_lineality; i++) {
00369       free(L[i]);
00370     }
00371     if (L != NULL) {
00372       free(L);
00373     }
00374   }
00375 
00376   for (i = 0; i < num_inequalities; i++) {
00377     free(Iq[i]);
00378   }
00379   free(Iq);
00380 }
00381 #endif
00382 
00383 #else
00384 
00385 int compute_cone_generators(int* n,...)
00386 {
00387   printf("Error: compute_cone_generators() should not be called without CDDLIB.\n\n");
00388   printf("\nAPPSPACK Error\n");
00389   exit(1);
00390 }
00391 
00392 #endif

 

© Sandia Corporation | Site Contact | Privacy and Security

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