Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (4 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/Clp_ampl.cpp

    r2280 r2385  
    3535
    3636#ifdef HAVE_UNISTD_H
    37 # include "unistd.h"
     37#include "unistd.h"
    3838#endif
    3939#include "CoinUtilsConfig.h"
     
    4646#include "Clp_ampl.h"
    4747extern "C" {
    48 # include "getstub.h"
    49 # include "asl_pfgh.h"
     48#include "getstub.h"
     49#include "asl_pfgh.h"
    5050}
    5151
     
    5353#include <cassert>
    5454/* so decodePhrase and clpCheck can access */
    55 static ampl_info * saveInfo = NULL;
     55static ampl_info *saveInfo = NULL;
    5656// Set to 1 if algorithm found
    5757static char algFound[20] = "";
    58 static char*
     58static char *
    5959checkPhrase(Option_Info *oi, keyword *kw, char *v)
    6060{
    61     if (strlen(v))
    62         printf("string %s\n", v);
    63     // Say algorithm found
    64     strcpy(algFound, kw->desc);
    65     return v;
    66 }
    67 static char*
     61  if (strlen(v))
     62    printf("string %s\n", v);
     63  // Say algorithm found
     64  strcpy(algFound, kw->desc);
     65  return v;
     66}
     67static char *
    6868checkPhrase2(Option_Info *oi, keyword *kw, char *v)
    6969{
    70     if (strlen(v))
    71         printf("string %s\n", v);
    72     // put out keyword
    73     saveInfo->arguments = (char **) realloc(saveInfo->arguments, (saveInfo->numberArguments + 1) * sizeof(char *));
    74     saveInfo->arguments[saveInfo->numberArguments++] = strdup(kw->desc);
    75     return v;
     70  if (strlen(v))
     71    printf("string %s\n", v);
     72  // put out keyword
     73  saveInfo->arguments = (char **)realloc(saveInfo->arguments, (saveInfo->numberArguments + 1) * sizeof(char *));
     74  saveInfo->arguments[saveInfo->numberArguments++] = strdup(kw->desc);
     75  return v;
    7676}
    7777static fint
    78 decodePhrase(char * phrase, ftnlen length)
    79 {
    80     char * blank = strchr(phrase, ' ');
    81     if (blank) {
    82         /* split arguments */
    83         *blank = '\0';
    84         saveInfo->arguments = (char **) realloc(saveInfo->arguments, (saveInfo->numberArguments + 2) * sizeof(char *));
    85         saveInfo->arguments[saveInfo->numberArguments++] = strdup(phrase);
    86         *blank = ' ';
    87         phrase = blank + 1; /* move on */
    88         if (strlen(phrase))
    89             saveInfo->arguments[saveInfo->numberArguments++] = strdup(phrase);
    90     } else if (strlen(phrase)) {
    91         saveInfo->arguments = (char **) realloc(saveInfo->arguments, (saveInfo->numberArguments + 1) * sizeof(char *));
    92         saveInfo->arguments[saveInfo->numberArguments++] = strdup(phrase);
    93     }
    94     return 0;
     78decodePhrase(char *phrase, ftnlen length)
     79{
     80  char *blank = strchr(phrase, ' ');
     81  if (blank) {
     82    /* split arguments */
     83    *blank = '\0';
     84    saveInfo->arguments = (char **)realloc(saveInfo->arguments, (saveInfo->numberArguments + 2) * sizeof(char *));
     85    saveInfo->arguments[saveInfo->numberArguments++] = strdup(phrase);
     86    *blank = ' ';
     87    phrase = blank + 1; /* move on */
     88    if (strlen(phrase))
     89      saveInfo->arguments[saveInfo->numberArguments++] = strdup(phrase);
     90  } else if (strlen(phrase)) {
     91    saveInfo->arguments = (char **)realloc(saveInfo->arguments, (saveInfo->numberArguments + 1) * sizeof(char *));
     92    saveInfo->arguments[saveInfo->numberArguments++] = strdup(phrase);
     93  }
     94  return 0;
    9595}
    9696static void
    97 sos_kludge(int nsos, int *sosbeg, double *sosref, int * sosind)
    98 {
    99     // Adjust sosref if necessary to make monotonic increasing
    100     int i, j, k;
    101     // first sort
    102     for (i = 0; i < nsos; i++) {
    103         k = sosbeg[i];
    104         int end = sosbeg[i+1];
    105         CoinSort_2(sosref + k, sosref + end, sosind + k);
    106     }
    107     double t, t1;
    108     for (i = j = 0; i++ < nsos; ) {
    109         k = sosbeg[i];
    110         t = sosref[j];
    111         while (++j < k) {
    112             t1 = sosref[j];
    113             t += 1e-10;
    114             if (t1 <= t)
    115                 sosref[j] = t1 = t + 1e-10;
    116             t = t1;
    117         }
    118     }
     97sos_kludge(int nsos, int *sosbeg, double *sosref, int *sosind)
     98{
     99  // Adjust sosref if necessary to make monotonic increasing
     100  int i, j, k;
     101  // first sort
     102  for (i = 0; i < nsos; i++) {
     103    k = sosbeg[i];
     104    int end = sosbeg[i + 1];
     105    CoinSort_2(sosref + k, sosref + end, sosind + k);
     106  }
     107  double t, t1;
     108  for (i = j = 0; i++ < nsos;) {
     109    k = sosbeg[i];
     110    t = sosref[j];
     111    while (++j < k) {
     112      t1 = sosref[j];
     113      t += 1e-10;
     114      if (t1 <= t)
     115        sosref[j] = t1 = t + 1e-10;
     116      t = t1;
     117    }
     118  }
    119119}
    120120static char xxxxxx[20];
    121 #define VP (char*)
     121#define VP (char *)
    122122static keyword keywds[] = { /* must be sorted */
    123     { const_cast<char*>("barrier"),  checkPhrase,  (char *) xxxxxx ,
    124         const_cast<char*>("-barrier")},
    125     { const_cast<char*>("dual"),     checkPhrase,  (char *) xxxxxx ,
    126       const_cast<char*>("-dualsimplex")},
    127     { const_cast<char*>("help"),     checkPhrase2, (char *) xxxxxx ,
    128       const_cast<char*>("-?")},
    129     { const_cast<char*>("initial"),  checkPhrase,  (char *) xxxxxx ,
    130       const_cast<char*>("-initialsolve")},
    131     { const_cast<char*>("max"),      checkPhrase2, (char *) xxxxxx ,
    132       const_cast<char*>("-maximize")},
    133     { const_cast<char*>("maximize"), checkPhrase2, (char *) xxxxxx ,
    134       const_cast<char*>("-maximize")},
    135     { const_cast<char*>("primal"),   checkPhrase,  (char *) xxxxxx ,
    136       const_cast<char*>("-primalsimplex")},
    137     { const_cast<char*>("quit"),     checkPhrase2, (char *) xxxxxx ,
    138       const_cast<char*>("-quit")},
    139     { const_cast<char*>("wantsol"),  WS_val,      NULL,
    140       const_cast<char*>("write .sol file (without -AMPL)")}
     123  { const_cast< char * >("barrier"), checkPhrase, (char *)xxxxxx,
     124    const_cast< char * >("-barrier") },
     125  { const_cast< char * >("dual"), checkPhrase, (char *)xxxxxx,
     126    const_cast< char * >("-dualsimplex") },
     127  { const_cast< char * >("help"), checkPhrase2, (char *)xxxxxx,
     128    const_cast< char * >("-?") },
     129  { const_cast< char * >("initial"), checkPhrase, (char *)xxxxxx,
     130    const_cast< char * >("-initialsolve") },
     131  { const_cast< char * >("max"), checkPhrase2, (char *)xxxxxx,
     132    const_cast< char * >("-maximize") },
     133  { const_cast< char * >("maximize"), checkPhrase2, (char *)xxxxxx,
     134    const_cast< char * >("-maximize") },
     135  { const_cast< char * >("primal"), checkPhrase, (char *)xxxxxx,
     136    const_cast< char * >("-primalsimplex") },
     137  { const_cast< char * >("quit"), checkPhrase2, (char *)xxxxxx,
     138    const_cast< char * >("-quit") },
     139  { const_cast< char * >("wantsol"), WS_val, NULL,
     140    const_cast< char * >("write .sol file (without -AMPL)") }
    141141};
    142142static Option_Info Oinfo = {
    143     const_cast<char*>("clp"),
    144     const_cast<char*>("CLP " CLP_VERSION),
    145     const_cast<char*>("clp_options"),
    146     keywds,
    147     nkeywds,
    148     0,
    149     0,
    150     0,
    151     decodePhrase,
    152     0,
    153     0,
    154     0,
    155     20130502
     143  const_cast< char * >("clp"),
     144  const_cast< char * >("CLP " CLP_VERSION),
     145  const_cast< char * >("clp_options"),
     146  keywds,
     147  nkeywds,
     148  0,
     149  0,
     150  0,
     151  decodePhrase,
     152  0,
     153  0,
     154  0,
     155  20130502
    156156};
    157157// strdup used to avoid g++ compiler warning
    158158static SufDecl suftab[] = {
    159159#ifdef JJF_ZERO
    160     { const_cast<char*>("current"), 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
    161     { const_cast<char*>("current"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
    162     { const_cast<char*>("direction"), 0, ASL_Sufkind_var },
    163     { const_cast<char*>("down"), 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
    164     { const_cast<char*>("down"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
    165     { const_cast<char*>("priority"), 0, ASL_Sufkind_var },
     160  { const_cast< char * >("current"), 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
     161  { const_cast< char * >("current"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
     162  { const_cast< char * >("direction"), 0, ASL_Sufkind_var },
     163  { const_cast< char * >("down"), 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
     164  { const_cast< char * >("down"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
     165  { const_cast< char * >("priority"), 0, ASL_Sufkind_var },
    166166#endif
    167     { const_cast<char*>("cut"), 0, ASL_Sufkind_con },
    168     { const_cast<char*>("direction"), 0, ASL_Sufkind_var },
    169     { const_cast<char*>("downPseudocost"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
    170     { const_cast<char*>("priority"), 0, ASL_Sufkind_var },
    171     { const_cast<char*>("ref"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
    172     { const_cast<char*>("sos"), 0, ASL_Sufkind_var },
    173     { const_cast<char*>("sos"), 0, ASL_Sufkind_con },
    174     { const_cast<char*>("sosno"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
    175     { const_cast<char*>("sosref"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
    176     { const_cast<char*>("special"), 0, ASL_Sufkind_var },
    177     { const_cast<char*>("special"), 0, ASL_Sufkind_con },
    178     /*{ const_cast<char*>("special"), 0, ASL_Sufkind_con },*/
    179     { const_cast<char*>("sstatus"), 0, ASL_Sufkind_var, 0 },
    180     { const_cast<char*>("sstatus"), 0, ASL_Sufkind_con, 0 },
    181     { const_cast<char*>("upPseudocost"), 0, ASL_Sufkind_var | ASL_Sufkind_real }
     167  { const_cast< char * >("cut"), 0, ASL_Sufkind_con },
     168  { const_cast< char * >("direction"), 0, ASL_Sufkind_var },
     169  { const_cast< char * >("downPseudocost"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
     170  { const_cast< char * >("priority"), 0, ASL_Sufkind_var },
     171  { const_cast< char * >("ref"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
     172  { const_cast< char * >("sos"), 0, ASL_Sufkind_var },
     173  { const_cast< char * >("sos"), 0, ASL_Sufkind_con },
     174  { const_cast< char * >("sosno"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
     175  { const_cast< char * >("sosref"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
     176  { const_cast< char * >("special"), 0, ASL_Sufkind_var },
     177  { const_cast< char * >("special"), 0, ASL_Sufkind_con },
     178  /*{ const_cast<char*>("special"), 0, ASL_Sufkind_con },*/
     179  { const_cast< char * >("sstatus"), 0, ASL_Sufkind_var, 0 },
     180  { const_cast< char * >("sstatus"), 0, ASL_Sufkind_con, 0 },
     181  { const_cast< char * >("upPseudocost"), 0, ASL_Sufkind_var | ASL_Sufkind_real }
    182182#ifdef JJF_ZERO
    183     { const_cast<char*>("unbdd"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly},
    184     { const_cast<char*>("up"), 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
    185     { const_cast<char*>("up"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly }
     183  { const_cast< char * >("unbdd"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
     184  { const_cast< char * >("up"), 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
     185  { const_cast< char * >("up"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly }
    186186#endif
    187187};
     
    194194mip_stuff(void)
    195195{
    196     int i;
    197     double *pseudoUp, *pseudoDown;
    198     int *priority, *direction;
    199     // To label cuts (there will be other uses for special)
    200     int *cut;
    201     // To label special variables - at present 1= must be >= 1 or <= -1
    202     int * special;
    203     SufDesc *dpup, *dpdown, *dpri, *ddir, *dcut, *dspecial;
    204 
    205     ddir = suf_get("direction", ASL_Sufkind_var);
    206     direction = ddir->u.i;
    207     dpri = suf_get("priority", ASL_Sufkind_var);
    208     priority = dpri->u.i;
    209     dspecial = suf_get("special", ASL_Sufkind_con);
    210     dcut = suf_get("cut", ASL_Sufkind_con);
     196  int i;
     197  double *pseudoUp, *pseudoDown;
     198  int *priority, *direction;
     199  // To label cuts (there will be other uses for special)
     200  int *cut;
     201  // To label special variables - at present 1= must be >= 1 or <= -1
     202  int *special;
     203  SufDesc *dpup, *dpdown, *dpri, *ddir, *dcut, *dspecial;
     204
     205  ddir = suf_get("direction", ASL_Sufkind_var);
     206  direction = ddir->u.i;
     207  dpri = suf_get("priority", ASL_Sufkind_var);
     208  priority = dpri->u.i;
     209  dspecial = suf_get("special", ASL_Sufkind_con);
     210  dcut = suf_get("cut", ASL_Sufkind_con);
     211  cut = dcut->u.i;
     212  if (!cut) {
     213    // try special
     214    dcut = suf_get("special", ASL_Sufkind_con);
    211215    cut = dcut->u.i;
    212     if (!cut) {
    213         // try special
    214         dcut = suf_get("special", ASL_Sufkind_con);
    215         cut = dcut->u.i;
    216     }
    217     dspecial = suf_get("special", ASL_Sufkind_var);
    218     special = dspecial->u.i;
    219     dpdown = suf_get("downPseudocost", ASL_Sufkind_var);
    220     pseudoDown = dpdown->u.r;
    221     dpup = suf_get("upPseudocost", ASL_Sufkind_var);
    222     pseudoUp = dpup->u.r;
    223     assert(saveInfo);
    224     int numberColumns = saveInfo->numberColumns;
    225     if (direction) {
    226         int baddir = 0;
    227         saveInfo->branchDirection = (int *) malloc(numberColumns * sizeof(int));
    228         for (i = 0; i < numberColumns; i++) {
    229             int value = direction[i];
    230             if (value < -1 || value > 1) {
    231                 baddir++;
    232                 value = 0;
    233             }
    234             saveInfo->branchDirection[i] = value;
     216  }
     217  dspecial = suf_get("special", ASL_Sufkind_var);
     218  special = dspecial->u.i;
     219  dpdown = suf_get("downPseudocost", ASL_Sufkind_var);
     220  pseudoDown = dpdown->u.r;
     221  dpup = suf_get("upPseudocost", ASL_Sufkind_var);
     222  pseudoUp = dpup->u.r;
     223  assert(saveInfo);
     224  int numberColumns = saveInfo->numberColumns;
     225  if (direction) {
     226    int baddir = 0;
     227    saveInfo->branchDirection = (int *)malloc(numberColumns * sizeof(int));
     228    for (i = 0; i < numberColumns; i++) {
     229      int value = direction[i];
     230      if (value < -1 || value > 1) {
     231        baddir++;
     232        value = 0;
     233      }
     234      saveInfo->branchDirection[i] = value;
     235    }
     236    if (baddir)
     237      fprintf(Stderr,
     238        "Treating %d .direction values outside [-1, 1] as 0.\n",
     239        baddir);
     240  }
     241  if (priority) {
     242    int badpri = 0;
     243    saveInfo->priorities = (int *)malloc(numberColumns * sizeof(int));
     244    for (i = 0; i < numberColumns; i++) {
     245      int value = priority[i];
     246      if (value < 0) {
     247        badpri++;
     248        value = 0;
     249      }
     250      saveInfo->priorities[i] = value;
     251    }
     252    if (badpri)
     253      fprintf(Stderr,
     254        "Treating %d negative .priority values as 0\n",
     255        badpri);
     256  }
     257  if (special) {
     258    int badspecial = 0;
     259    saveInfo->special = (int *)malloc(numberColumns * sizeof(int));
     260    for (i = 0; i < numberColumns; i++) {
     261      int value = special[i];
     262      if (value < 0) {
     263        badspecial++;
     264        value = 0;
     265      }
     266      saveInfo->special[i] = value;
     267    }
     268    if (badspecial)
     269      fprintf(Stderr,
     270        "Treating %d negative special values as 0\n",
     271        badspecial);
     272  }
     273  int numberRows = saveInfo->numberRows;
     274  if (cut) {
     275    int badcut = 0;
     276    saveInfo->cut = (int *)malloc(numberRows * sizeof(int));
     277    for (i = 0; i < numberRows; i++) {
     278      int value = cut[i];
     279      if (value < 0) {
     280        badcut++;
     281        value = 0;
     282      }
     283      saveInfo->cut[i] = value;
     284    }
     285    if (badcut)
     286      fprintf(Stderr,
     287        "Treating %d negative cut values as 0\n",
     288        badcut);
     289  }
     290  if (pseudoDown || pseudoUp) {
     291    int badpseudo = 0;
     292    if (!pseudoDown || !pseudoUp)
     293      fprintf(Stderr,
     294        "Only one set of pseudocosts - assumed same\n");
     295    saveInfo->pseudoDown = (double *)malloc(numberColumns * sizeof(double));
     296    saveInfo->pseudoUp = (double *)malloc(numberColumns * sizeof(double));
     297    for (i = 0; i < numberColumns; i++) {
     298      double valueD = 0.0, valueU = 0.0;
     299      if (pseudoDown) {
     300        valueD = pseudoDown[i];
     301        if (valueD < 0) {
     302          badpseudo++;
     303          valueD = 0.0;
    235304        }
    236         if (baddir)
    237             fprintf(Stderr,
    238                     "Treating %d .direction values outside [-1, 1] as 0.\n",
    239                     baddir);
    240     }
    241     if (priority) {
    242         int badpri = 0;
    243         saveInfo->priorities = (int *) malloc(numberColumns * sizeof(int));
    244         for (i = 0; i < numberColumns; i++) {
    245             int value = priority[i];
    246             if (value < 0) {
    247                 badpri++;
    248                 value = 0;
    249             }
    250             saveInfo->priorities[i] = value;
     305      }
     306      if (pseudoUp) {
     307        valueU = pseudoUp[i];
     308        if (valueU < 0) {
     309          badpseudo++;
     310          valueU = 0.0;
    251311        }
    252         if (badpri)
    253             fprintf(Stderr,
    254                     "Treating %d negative .priority values as 0\n",
    255                     badpri);
    256     }
    257     if (special) {
    258         int badspecial = 0;
    259         saveInfo->special = (int *) malloc(numberColumns * sizeof(int));
    260         for (i = 0; i < numberColumns; i++) {
    261             int value = special[i];
    262             if (value < 0) {
    263                 badspecial++;
    264                 value = 0;
    265             }
    266             saveInfo->special[i] = value;
    267         }
    268         if (badspecial)
    269             fprintf(Stderr,
    270                     "Treating %d negative special values as 0\n",
    271                     badspecial);
    272     }
    273     int numberRows = saveInfo->numberRows;
    274     if (cut) {
    275         int badcut = 0;
    276         saveInfo->cut = (int *) malloc(numberRows * sizeof(int));
    277         for (i = 0; i < numberRows; i++) {
    278             int value = cut[i];
    279             if (value < 0) {
    280                 badcut++;
    281                 value = 0;
    282             }
    283             saveInfo->cut[i] = value;
    284         }
    285         if (badcut)
    286             fprintf(Stderr,
    287                     "Treating %d negative cut values as 0\n",
    288                     badcut);
    289     }
    290     if (pseudoDown || pseudoUp) {
    291         int badpseudo = 0;
    292         if (!pseudoDown || !pseudoUp)
    293             fprintf(Stderr,
    294                     "Only one set of pseudocosts - assumed same\n");
    295         saveInfo->pseudoDown = (double *) malloc(numberColumns * sizeof(double));
    296         saveInfo->pseudoUp = (double *) malloc(numberColumns * sizeof(double));
    297         for (i = 0; i < numberColumns; i++) {
    298             double valueD = 0.0, valueU = 0.0;
    299             if (pseudoDown) {
    300                 valueD = pseudoDown[i];
    301                 if (valueD < 0) {
    302                     badpseudo++;
    303                     valueD = 0.0;
    304                 }
    305             }
    306             if (pseudoUp) {
    307                 valueU = pseudoUp[i];
    308                 if (valueU < 0) {
    309                     badpseudo++;
    310                     valueU = 0.0;
    311                 }
    312             }
    313             if (!valueD)
    314                 valueD = valueU;
    315             if (!valueU)
    316                 valueU = valueD;
    317             saveInfo->pseudoDown[i] = valueD;
    318             saveInfo->pseudoUp[i] = valueU;
    319         }
    320         if (badpseudo)
    321             fprintf(Stderr,
    322                     "Treating %d negative pseudoCosts as 0.0\n", badpseudo);
    323     }
     312      }
     313      if (!valueD)
     314        valueD = valueU;
     315      if (!valueU)
     316        valueU = valueD;
     317      saveInfo->pseudoDown[i] = valueD;
     318      saveInfo->pseudoUp[i] = valueU;
     319    }
     320    if (badpseudo)
     321      fprintf(Stderr,
     322        "Treating %d negative pseudoCosts as 0.0\n", badpseudo);
     323  }
    324324}
    325325static void
    326326stat_map(int *stat, int n, int *map, int mx, const char *what)
    327327{
    328     int bad, i, i1 = 0, j, j1 = 0;
    329     static char badfmt[] = "Coin driver: %s[%d] = %d\n";
    330 
    331     for (i = bad = 0; i < n; i++) {
    332         if ((j = stat[i]) >= 0 && j <= mx)
    333             stat[i] = map[j];
    334         else {
    335             stat[i] = 0;
    336             i1 = i;
    337             j1 = j;
    338             if (!bad++)
    339                 fprintf(Stderr, badfmt, what, i, j);
     328  int bad, i, i1 = 0, j, j1 = 0;
     329  static char badfmt[] = "Coin driver: %s[%d] = %d\n";
     330
     331  for (i = bad = 0; i < n; i++) {
     332    if ((j = stat[i]) >= 0 && j <= mx)
     333      stat[i] = map[j];
     334    else {
     335      stat[i] = 0;
     336      i1 = i;
     337      j1 = j;
     338      if (!bad++)
     339        fprintf(Stderr, badfmt, what, i, j);
     340    }
     341  }
     342  if (bad > 1) {
     343    if (bad == 2)
     344      fprintf(Stderr, badfmt, what, i1, j1);
     345    else
     346      fprintf(Stderr,
     347        "Coin driver: %d messages about bad %s values suppressed.\n",
     348        bad - 1, what);
     349  }
     350}
     351
     352int readAmpl(ampl_info *info, int argc, char **argv, void **coinModel)
     353{
     354  char *stub;
     355  ograd *og;
     356  int i;
     357  SufDesc *csd;
     358  SufDesc *rsd;
     359  /*bool *basis, *lower;*/
     360  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
     361  char *environment = getenv("clp_options");
     362  char tempBuffer[20];
     363  double *obj;
     364  double *columnLower;
     365  double *columnUpper;
     366  double *rowLower;
     367  double *rowUpper;
     368  char **saveArgv = argv;
     369  char fileName[1000];
     370  if (argc > 1)
     371    strcpy(fileName, argv[1]);
     372  else
     373    fileName[0] = '\0';
     374  int nonLinearType = -1;
     375  // testosi parameter - if >= 10 then go in through coinModel
     376  for (i = 1; i < argc; i++) {
     377    if (!strncmp(argv[i], "testosi", 7)) {
     378      char *equals = strchr(argv[i], '=');
     379      if (equals && atoi(equals + 1) >= 10 && atoi(equals + 1) <= 20) {
     380        nonLinearType = atoi(equals + 1);
     381        break;
     382      }
     383    }
     384  }
     385  int saveArgc = argc;
     386  if (info->numberRows != -1234567)
     387    memset(info, 0, sizeof(ampl_info)); // overwrite unless magic number set
     388  /* save so can be accessed by decodePhrase */
     389  saveInfo = info;
     390  info->numberArguments = 0;
     391  info->arguments = (char **)malloc(2 * sizeof(char *));
     392  info->arguments[info->numberArguments++] = strdup("ampl");
     393  info->arguments[info->numberArguments++] = strdup("clp");
     394  asl = ASL_alloc(ASL_read_f);
     395  stub = getstub(&argv, &Oinfo);
     396  if (!stub)
     397    usage_ASL(&Oinfo, 1);
     398  nl = jac0dim(stub, 0);
     399  suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
     400
     401  /* set A_vals to get the constraints column-wise (malloc so can be freed) */
     402  A_vals = (double *)malloc(nzc * sizeof(double));
     403  if (!A_vals) {
     404    printf("no memory\n");
     405    return 1;
     406  }
     407  /* say we want primal solution */
     408  want_xpi0 = 1;
     409  /* for basis info */
     410  info->columnStatus = (int *)malloc(n_var * sizeof(int));
     411  for (int i = 0; i < n_var; i++)
     412    info->columnStatus[i] = 3;
     413  info->rowStatus = (int *)malloc(n_con * sizeof(int));
     414  for (int i = 0; i < n_con; i++)
     415    info->rowStatus[i] = 1;
     416  csd = suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
     417  rsd = suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
     418  if (!(nlvc + nlvo) && nonLinearType < 10) {
     419    /* read linear model*/
     420#if COIN_BIG_INDEX == 2
     421    f_read(nl, ASL_allow_Z | ASL_use_Z);
     422#else
     423    f_read(nl, 0);
     424#endif
     425    // see if any sos
     426    if (true) {
     427      char *sostype;
     428      int nsosnz, *sosbeg, *sosind, *sospri;
     429      double *sosref;
     430      int nsos;
     431      int i = ASL_suf_sos_explict_free;
     432      int copri[2], **p_sospri;
     433      copri[0] = 0;
     434      copri[1] = 0;
     435      p_sospri = &sospri;
     436      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
     437        &sosbeg, &sosind, &sosref);
     438      if (nsos) {
     439        info->numberSos = nsos;
     440        info->sosType = (char *)malloc(nsos);
     441        info->sosPriority = (int *)malloc(nsos * sizeof(int));
     442        info->sosStart = (int *)malloc((nsos + 1) * sizeof(int));
     443        info->sosIndices = (int *)malloc(nsosnz * sizeof(int));
     444        info->sosReference = (double *)malloc(nsosnz * sizeof(double));
     445        sos_kludge(nsos, sosbeg, sosref, sosind);
     446        for (int i = 0; i < nsos; i++) {
     447          char ichar = sostype[i];
     448          assert(ichar == '1' || ichar == '2');
     449          info->sosType[i] = static_cast< char >(ichar - '0');
    340450        }
    341     }
    342     if (bad > 1) {
    343         if (bad == 2)
    344             fprintf(Stderr, badfmt, what, i1, j1);
    345         else
    346             fprintf(Stderr,
    347                     "Coin driver: %d messages about bad %s values suppressed.\n",
    348                     bad - 1, what);
    349     }
    350 }
    351 
    352 int
    353 readAmpl(ampl_info * info, int argc, char **argv, void ** coinModel)
    354 {
    355     char *stub;
    356     ograd *og;
    357     int i;
    358     SufDesc *csd;
    359     SufDesc *rsd;
    360     /*bool *basis, *lower;*/
    361     /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
    362     char * environment = getenv("clp_options");
    363     char tempBuffer[20];
    364     double * obj;
    365     double * columnLower;
    366     double * columnUpper;
    367     double * rowLower;
    368     double * rowUpper;
    369     char ** saveArgv = argv;
    370     char fileName[1000];
    371     if (argc > 1)
    372         strcpy(fileName, argv[1]);
     451        memcpy(info->sosPriority, sospri, nsos * sizeof(int));
     452        memcpy(info->sosStart, sosbeg, (nsos + 1) * sizeof(int));
     453        memcpy(info->sosIndices, sosind, nsosnz * sizeof(int));
     454        memcpy(info->sosReference, sosref, nsosnz * sizeof(double));
     455      }
     456    }
     457
     458    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
     459    Oinfo.uinfo = tempBuffer;
     460    if (getopts(argv, &Oinfo))
     461      return 1;
     462    /* objective*/
     463    obj = (double *)malloc(n_var * sizeof(double));
     464    for (i = 0; i < n_var; i++)
     465      obj[i] = 0.0;
     466    if (n_obj) {
     467      for (og = Ograd[0]; og; og = og->next)
     468        obj[og->varno] = og->coef;
     469    }
     470    if (objtype[0])
     471      info->direction = -1.0;
    373472    else
    374         fileName[0] = '\0';
    375     int nonLinearType = -1;
    376     // testosi parameter - if >= 10 then go in through coinModel
    377     for (i = 1; i < argc; i++) {
    378         if (!strncmp(argv[i], "testosi", 7)) {
    379             char * equals = strchr(argv[i], '=');
    380             if (equals && atoi(equals + 1) >= 10 && atoi(equals + 1) <= 20) {
    381                 nonLinearType = atoi(equals + 1);
    382                 break;
    383             }
     473      info->direction = 1.0;
     474    info->offset = objconst(0);
     475    /* Column bounds*/
     476    columnLower = (double *)malloc(n_var * sizeof(double));
     477    columnUpper = (double *)malloc(n_var * sizeof(double));
     478    for (i = 0; i < n_var; i++) {
     479      columnLower[i] = LUv[2 * i];
     480      if (columnLower[i] <= negInfinity)
     481        columnLower[i] = -COIN_DBL_MAX;
     482      columnUpper[i] = LUv[2 * i + 1];
     483      if (columnUpper[i] >= Infinity)
     484        columnUpper[i] = COIN_DBL_MAX;
     485    }
     486    /* Row bounds*/
     487    rowLower = (double *)malloc(n_con * sizeof(double));
     488    rowUpper = (double *)malloc(n_con * sizeof(double));
     489    for (i = 0; i < n_con; i++) {
     490      rowLower[i] = LUrhs[2 * i];
     491      if (rowLower[i] <= negInfinity)
     492        rowLower[i] = -COIN_DBL_MAX;
     493      rowUpper[i] = LUrhs[2 * i + 1];
     494      if (rowUpper[i] >= Infinity)
     495        rowUpper[i] = COIN_DBL_MAX;
     496    }
     497    info->numberRows = n_con;
     498    info->numberColumns = n_var;
     499    info->numberElements = nzc;
     500    info->numberBinary = nbv;
     501    info->numberIntegers = niv + nbv;
     502    info->objective = obj;
     503    info->rowLower = rowLower;
     504    info->rowUpper = rowUpper;
     505    info->columnLower = columnLower;
     506    info->columnUpper = columnUpper;
     507#if COIN_BIG_INDEX == 0
     508    info->starts = A_colstarts;
     509#else
     510    info->starts = A_colstartsZ;
     511#endif
     512    /*A_colstarts=NULL;*/
     513    info->rows = A_rownos;
     514    /*A_rownos=NULL;*/
     515    info->elements = A_vals;
     516    /*A_vals=NULL;*/
     517    info->primalSolution = NULL;
     518    /* put in primalSolution if exists */
     519    if (X0) {
     520      info->primalSolution = (double *)malloc(n_var * sizeof(double));
     521      memcpy(info->primalSolution, X0, n_var * sizeof(double));
     522    }
     523    info->dualSolution = NULL;
     524    if (niv + nbv > 0)
     525      mip_stuff(); // get any extra info
     526    if ((!(niv + nbv) && (csd->kind & ASL_Sufkind_input))
     527      || (rsd->kind & ASL_Sufkind_input)) {
     528      /* convert status - need info on map */
     529      static int map[] = { 1, 3, 1, 1, 2, 1, 1 };
     530      stat_map(info->columnStatus, n_var, map, 6, "incoming columnStatus");
     531      stat_map(info->rowStatus, n_con, map, 6, "incoming rowStatus");
     532    } else {
     533      /* all slack basis */
     534      // leave status for output */
     535#ifdef JJF_ZERO
     536      free(info->rowStatus);
     537      info->rowStatus = NULL;
     538      free(info->columnStatus);
     539      info->columnStatus = NULL;
     540#endif
     541    }
     542  } else {
     543    // QP
     544    // Add .nl if not there
     545    if (!strstr(fileName, ".nl"))
     546      strcat(fileName, ".nl");
     547    CoinModel *model = new CoinModel((nonLinearType > 10) ? 2 : 1, fileName, info);
     548    if (model->numberRows() > 0 || model->numberColumns() > 0)
     549      *coinModel = (void *)model;
     550    Oinfo.uinfo = tempBuffer;
     551    if (getopts(argv, &Oinfo))
     552      return 1;
     553    Oinfo.wantsol = 1;
     554    if (objtype[0])
     555      info->direction = -1.0;
     556    else
     557      info->direction = 1.0;
     558    model->setOptimizationDirection(info->direction);
     559    info->offset = objconst(0);
     560    info->numberRows = n_con;
     561    info->numberColumns = n_var;
     562    info->numberElements = nzc;
     563    info->numberBinary = nbv;
     564    int numberIntegers = niv + nlvci + nlvoi + nbv;
     565    if (nlvci + nlvoi + nlvc + nlvo) {
     566      // Non linear
     567      // No idea if there are overlaps so compute
     568      int numberIntegers = 0;
     569      for (i = 0; i < n_var; i++) {
     570        if (model->columnIsInteger(i))
     571          numberIntegers++;
     572      }
     573    }
     574    info->numberIntegers = numberIntegers;
     575    // Say nonlinear if it is
     576    info->nonLinear = nlvc + nlvo;
     577    if (numberIntegers > 0) {
     578      mip_stuff(); // get any extra info
     579      if (info->cut)
     580        model->setCutMarker(info->numberRows, info->cut);
     581      if (info->priorities)
     582        model->setPriorities(info->numberColumns, info->priorities);
     583    }
     584  }
     585  /* add -solve - unless something there already
     586     - also check for sleep=yes */
     587  {
     588    int found = 0;
     589    int foundLog = 0;
     590    int foundSleep = 0;
     591    const char *something[] = { "solve", "branch", "duals", "primals", "user" };
     592    for (i = 0; i < info->numberArguments; i++) {
     593      unsigned int j;
     594      const char *argument = info->arguments[i];
     595      for (j = 0; j < sizeof(something) / sizeof(char *); j++) {
     596        const char *check = something[j];
     597        if (!strncmp(argument, check, sizeof(check))) {
     598          found = (int)(j + 1);
     599        } else if (!strncmp(argument, "log", 3)) {
     600          foundLog = 1;
     601        } else if (!strncmp(argument, "sleep", 5)) {
     602          foundSleep = 1;
    384603        }
    385     }
    386     int saveArgc = argc;
    387     if (info->numberRows != -1234567)
    388         memset(info, 0, sizeof(ampl_info)); // overwrite unless magic number set
    389     /* save so can be accessed by decodePhrase */
    390     saveInfo = info;
    391     info->numberArguments = 0;
    392     info->arguments = (char **) malloc(2 * sizeof(char *));
    393     info->arguments[info->numberArguments++] = strdup("ampl");
    394     info->arguments[info->numberArguments++] = strdup("clp");
     604      }
     605    }
     606    if (foundLog) {
     607      /* print options etc */
     608      for (i = 0; i < saveArgc; i++)
     609        printf("%s ", saveArgv[i]);
     610      printf("\n");
     611      if (environment)
     612        printf("env %s\n", environment);
     613      /*printf("%d rows %d columns %d elements\n",n_con,n_var,nzc);*/
     614    }
     615    if (!found) {
     616      if (!strlen(algFound)) {
     617        info->arguments = (char **)realloc(info->arguments, (info->numberArguments + 1) * sizeof(char *));
     618        info->arguments[info->numberArguments++] = strdup("-solve");
     619      } else {
     620        // use algorithm from keyword
     621        info->arguments = (char **)realloc(info->arguments, (info->numberArguments + 1) * sizeof(char *));
     622        info->arguments[info->numberArguments++] = strdup(algFound);
     623      }
     624    }
     625    if (foundSleep) {
     626      /* let user copy .nl file */
     627      fprintf(stderr, "You can copy .nl file %s for debug purposes or attach debugger\n", saveArgv[1]);
     628      fprintf(stderr, "Type q to quit, anything else to continue\n");
     629      int getChar = getc(stdin);
     630      if (getChar == 'q' || getChar == 'Q')
     631        exit(1);
     632    }
     633  }
     634  /* add -quit */
     635  info->arguments = (char **)realloc(info->arguments, (info->numberArguments + 1) * sizeof(char *));
     636  info->arguments[info->numberArguments++] = strdup("-quit");
     637  return 0;
     638}
     639void freeArrays1(ampl_info *info)
     640{
     641  free(info->objective);
     642  info->objective = NULL;
     643  free(info->rowLower);
     644  info->rowLower = NULL;
     645  free(info->rowUpper);
     646  info->rowUpper = NULL;
     647  free(info->columnLower);
     648  info->columnLower = NULL;
     649  free(info->columnUpper);
     650  info->columnUpper = NULL;
     651  /* this one not freed by ASL_free */
     652  free(info->elements);
     653  info->elements = NULL;
     654  free(info->primalSolution);
     655  info->primalSolution = NULL;
     656  free(info->dualSolution);
     657  info->dualSolution = NULL;
     658  /*free(info->rowStatus);
     659    info->rowStatus=NULL;
     660    free(info->columnStatus);
     661    info->columnStatus=NULL;*/
     662}
     663void freeArrays2(ampl_info *info)
     664{
     665  free(info->primalSolution);
     666  info->primalSolution = NULL;
     667  free(info->dualSolution);
     668  info->dualSolution = NULL;
     669  free(info->rowStatus);
     670  info->rowStatus = NULL;
     671  free(info->columnStatus);
     672  info->columnStatus = NULL;
     673  free(info->priorities);
     674  info->priorities = NULL;
     675  free(info->branchDirection);
     676  info->branchDirection = NULL;
     677  free(info->pseudoDown);
     678  info->pseudoDown = NULL;
     679  free(info->pseudoUp);
     680  info->pseudoUp = NULL;
     681  free(info->sosType);
     682  info->sosType = NULL;
     683  free(info->sosPriority);
     684  info->sosPriority = NULL;
     685  free(info->sosStart);
     686  info->sosStart = NULL;
     687  free(info->sosIndices);
     688  info->sosIndices = NULL;
     689  free(info->sosReference);
     690  info->sosReference = NULL;
     691  free(info->cut);
     692  info->cut = NULL;
     693  ASL_free(&asl);
     694}
     695void freeArgs(ampl_info *info)
     696{
     697  int i;
     698  for (i = 0; i < info->numberArguments; i++)
     699    free(info->arguments[i]);
     700  free(info->arguments);
     701}
     702int ampl_obj_prec()
     703{
     704  int precision = obj_prec();
     705  if (precision <= 0)
     706    precision = 15;
     707  return precision;
     708}
     709void writeAmpl(ampl_info *info)
     710{
     711  char buf[1000];
     712  typedef struct {
     713    const char *msg;
     714    int code;
     715    int wantObj;
     716  } Sol_info;
     717  static Sol_info solinfo[] = {
     718    { "optimal solution", 000, 1 },
     719    { "infeasible", 200, 1 },
     720    { "unbounded", 300, 0 },
     721    { "iteration limit etc", 400, 1 },
     722    { "solution limit", 401, 1 },
     723    { "ran out of space", 500, 0 },
     724    { "status unknown", 501, 1 },
     725    { "bug!", 502, 0 },
     726    { "best MIP solution so far restored", 101, 1 },
     727    { "failed to restore best MIP solution", 503, 1 },
     728    { "optimal (?) solution", 100, 1 }
     729  };
     730  /* convert status - need info on map */
     731  static int map[] = { 0, 3, 4, 1 };
     732  sprintf(buf, "%s %s", Oinfo.bsname, info->buffer);
     733  solve_result_num = solinfo[info->problemStatus].code;
     734  if (info->columnStatus) {
     735    stat_map(info->columnStatus, n_var, map, 4, "outgoing columnStatus");
     736    stat_map(info->rowStatus, n_con, map, 4, "outgoing rowStatus");
     737    suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
     738    suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
     739  }
     740  write_sol(buf, info->primalSolution, info->dualSolution, &Oinfo);
     741}
     742/* Read a problem from AMPL nl file
     743 */
     744CoinModel::CoinModel(int nonLinear, const char *fileName, const void *info)
     745  : CoinBaseModel()
     746  , maximumRows_(0)
     747  , maximumColumns_(0)
     748  , numberElements_(0)
     749  , maximumElements_(0)
     750  , numberQuadraticElements_(0)
     751  , maximumQuadraticElements_(0)
     752  , rowLower_(NULL)
     753  , rowUpper_(NULL)
     754  , rowType_(NULL)
     755  , objective_(NULL)
     756  , columnLower_(NULL)
     757  , columnUpper_(NULL)
     758  , integerType_(NULL)
     759  , columnType_(NULL)
     760  , start_(NULL)
     761  , elements_(NULL)
     762  , packedMatrix_(NULL)
     763  , quadraticElements_(NULL)
     764  , sortIndices_(NULL)
     765  , sortElements_(NULL)
     766  , sortSize_(0)
     767  , sizeAssociated_(0)
     768  , associated_(NULL)
     769  , numberSOS_(0)
     770  , startSOS_(NULL)
     771  , memberSOS_(NULL)
     772  , typeSOS_(NULL)
     773  , prioritySOS_(NULL)
     774  , referenceSOS_(NULL)
     775  , priority_(NULL)
     776  , cut_(NULL)
     777  , moreInfo_(NULL)
     778  , type_(-1)
     779  , noNames_(false)
     780  , links_(0)
     781{
     782  problemName_ = "";
     783  int status = 0;
     784  if (!strcmp(fileName, "-") || !strcmp(fileName, "stdin")) {
     785    // stdin
     786  } else {
     787    std::string name = fileName;
     788    bool readable = fileCoinReadable(name);
     789    if (!readable) {
     790      std::cerr << "Unable to open file "
     791                << fileName << std::endl;
     792      status = -1;
     793    }
     794  }
     795  if (!status) {
     796    gdb(nonLinear, fileName, info);
     797  }
     798}
     799#ifdef JJF_ZERO
     800static real
     801qterm(ASL *asl, fint *colq, fint *rowq, real *delsq)
     802{
     803  double t, t1, *x, *x0, *xe;
     804  fint *rq0, *rqe;
     805
     806  t = 0.;
     807  x0 = x = X0;
     808  xe = x + n_var;
     809  rq0 = rowq;
     810  while (x < xe) {
     811    t1 = *x++;
     812    rqe = rq0 + *++colq;
     813    while (rowq < rqe)
     814      t += t1 * x0[*rowq++] * *delsq++;
     815  }
     816  return 0.5 * t;
     817}
     818#endif
     819// stolen from IPopt with changes
     820typedef struct {
     821  double obj_sign_;
     822  ASL_pfgh *asl_;
     823  double *non_const_x_;
     824  int *column_; // for jacobian
     825  int *rowStart_;
     826  double *gradient_;
     827  double *constraintValues_;
     828  int nz_h_full_; // number of nonzeros in hessian
     829  int nerror_;
     830  bool objval_called_with_current_x_;
     831  bool conval_called_with_current_x_;
     832  bool jacval_called_with_current_x_;
     833} ClpAmplInfo;
     834
     835void CoinModel::gdb(int nonLinear, const char *fileName, const void *info)
     836{
     837  const ampl_info *amplInfo = (const ampl_info *)info;
     838  ograd *og = NULL;
     839  int i;
     840  SufDesc *csd = NULL;
     841  SufDesc *rsd = NULL;
     842  /*bool *basis, *lower;*/
     843  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
     844  //char tempBuffer[20];
     845  double *objective = NULL;
     846  double *columnLower = NULL;
     847  double *columnUpper = NULL;
     848  double *rowLower = NULL;
     849  double *rowUpper = NULL;
     850  int *columnStatus = NULL;
     851  int *rowStatus = NULL;
     852  int numberRows = -1;
     853  int numberColumns = -1;
     854  int numberElements = -1;
     855  int numberBinary = -1;
     856  int numberIntegers = -1;
     857  int numberAllNonLinearBoth = 0;
     858  int numberIntegerNonLinearBoth = 0;
     859  int numberAllNonLinearConstraints = 0;
     860  int numberIntegerNonLinearConstraints = 0;
     861  int numberAllNonLinearObjective = 0;
     862  int numberIntegerNonLinearObjective = 0;
     863  double *primalSolution = NULL;
     864  double direction = 1.0;
     865  char *stub = strdup(fileName);
     866  CoinPackedMatrix matrixByRow;
     867  fint **colqp = NULL;
     868  int *z = NULL;
     869  if (nonLinear == 0) {
     870    // linear
    395871    asl = ASL_alloc(ASL_read_f);
    396     stub = getstub(&argv, &Oinfo);
    397     if (!stub)
    398         usage_ASL(&Oinfo, 1);
    399872    nl = jac0dim(stub, 0);
     873    free(stub);
    400874    suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
    401875
    402876    /* set A_vals to get the constraints column-wise (malloc so can be freed) */
    403     A_vals = (double *) malloc(nzc * sizeof(double));
     877    A_vals = (double *)malloc(nzc * sizeof(double));
    404878    if (!A_vals) {
    405         printf("no memory\n");
    406         return 1;
     879      printf("no memory\n");
     880      return;
    407881    }
    408882    /* say we want primal solution */
    409883    want_xpi0 = 1;
    410884    /* for basis info */
    411     info->columnStatus = (int *) malloc(n_var * sizeof(int));
    412     for (int i=0;i<n_var;i++)
    413       info->columnStatus[i]=3;
    414     info->rowStatus = (int *) malloc(n_con * sizeof(int));
    415     for (int i=0;i<n_con;i++)
    416       info->rowStatus[i]=1;
    417     csd = suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
    418     rsd = suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
    419     if (!(nlvc + nlvo) && nonLinearType < 10) {
    420         /* read linear model*/
    421 #if COIN_BIG_INDEX==2
    422         f_read(nl, ASL_allow_Z|ASL_use_Z);
     885    columnStatus = (int *)malloc(n_var * sizeof(int));
     886    rowStatus = (int *)malloc(n_con * sizeof(int));
     887    csd = suf_iput("sstatus", ASL_Sufkind_var, columnStatus);
     888    rsd = suf_iput("sstatus", ASL_Sufkind_con, rowStatus);
     889    /* read linear model*/
     890#if COIN_BIG_INDEX == 2
     891    f_read(nl, ASL_allow_Z | ASL_use_Z);
    423892#else
    424         f_read(nl, 0);
     893    f_read(nl, 0);
    425894#endif
    426         // see if any sos
    427         if (true) {
    428             char *sostype;
    429             int nsosnz, *sosbeg, *sosind, * sospri;
    430             double *sosref;
    431             int nsos;
    432             int i = ASL_suf_sos_explict_free;
    433             int copri[2], **p_sospri;
    434             copri[0] = 0;
    435             copri[1] = 0;
    436             p_sospri = &sospri;
    437             nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
    438                            &sosbeg, &sosind, &sosref);
    439             if (nsos) {
    440                 info->numberSos = nsos;
    441                 info->sosType = (char *) malloc(nsos);
    442                 info->sosPriority = (int *) malloc(nsos * sizeof(int));
    443                 info->sosStart = (int *) malloc((nsos + 1) * sizeof(int));
    444                 info->sosIndices = (int *) malloc(nsosnz * sizeof(int));
    445                 info->sosReference = (double *) malloc(nsosnz * sizeof(double));
    446                 sos_kludge(nsos, sosbeg, sosref, sosind);
    447                 for (int i = 0; i < nsos; i++) {
    448                     char ichar = sostype[i];
    449                     assert (ichar == '1' || ichar == '2');
    450                     info->sosType[i] = static_cast<char>(ichar - '0');
    451                 }
    452                 memcpy(info->sosPriority, sospri, nsos*sizeof(int));
    453                 memcpy(info->sosStart, sosbeg, (nsos + 1)*sizeof(int));
    454                 memcpy(info->sosIndices, sosind, nsosnz*sizeof(int));
    455                 memcpy(info->sosReference, sosref, nsosnz*sizeof(double));
     895    // see if any sos
     896    if (true) {
     897      char *sostype;
     898      int nsosnz, *sosbeg, *sosind, *sospri;
     899      double *sosref;
     900      int nsos;
     901      int i = ASL_suf_sos_explict_free;
     902      int copri[2], **p_sospri;
     903      copri[0] = 0;
     904      copri[1] = 0;
     905      p_sospri = &sospri;
     906      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
     907        &sosbeg, &sosind, &sosref);
     908      if (nsos) {
     909        abort();
     910#ifdef JJF_ZERO
     911        info->numberSos = nsos;
     912        info->sosType = (char *)malloc(nsos);
     913        info->sosPriority = (int *)malloc(nsos * sizeof(int));
     914        info->sosStart = (int *)malloc((nsos + 1) * sizeof(int));
     915        info->sosIndices = (int *)malloc(nsosnz * sizeof(int));
     916        info->sosReference = (double *)malloc(nsosnz * sizeof(double));
     917        sos_kludge(nsos, sosbeg, sosref, sosind);
     918        for (int i = 0; i < nsos; i++) {
     919          int ichar = sostype[i];
     920          assert(ichar == '1' || ichar == '2');
     921          info->sosType[i] = ichar - '0';
     922        }
     923        memcpy(info->sosPriority, sospri, nsos * sizeof(int));
     924        memcpy(info->sosStart, sosbeg, (nsos + 1) * sizeof(int));
     925        memcpy(info->sosIndices, sosind, nsosnz * sizeof(int));
     926        memcpy(info->sosReference, sosref, nsosnz * sizeof(double));
     927#endif
     928      }
     929    }
     930
     931    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
     932    //Oinfo.uinfo = tempBuffer;
     933    //if (getopts(argv, &Oinfo))
     934    //return 1;
     935    /* objective*/
     936    objective = (double *)malloc(n_var * sizeof(double));
     937    for (i = 0; i < n_var; i++)
     938      objective[i] = 0.0;
     939    if (n_obj) {
     940      for (og = Ograd[0]; og; og = og->next)
     941        objective[og->varno] = og->coef;
     942    }
     943    if (objtype[0])
     944      direction = -1.0;
     945    else
     946      direction = 1.0;
     947    objectiveOffset_ = objconst(0);
     948    /* Column bounds*/
     949    columnLower = (double *)malloc(n_var * sizeof(double));
     950    columnUpper = (double *)malloc(n_var * sizeof(double));
     951    for (i = 0; i < n_var; i++) {
     952      columnLower[i] = LUv[2 * i];
     953      if (columnLower[i] <= negInfinity)
     954        columnLower[i] = -COIN_DBL_MAX;
     955      columnUpper[i] = LUv[2 * i + 1];
     956      if (columnUpper[i] >= Infinity)
     957        columnUpper[i] = COIN_DBL_MAX;
     958    }
     959    /* Row bounds*/
     960    rowLower = (double *)malloc(n_con * sizeof(double));
     961    rowUpper = (double *)malloc(n_con * sizeof(double));
     962    for (i = 0; i < n_con; i++) {
     963      rowLower[i] = LUrhs[2 * i];
     964      if (rowLower[i] <= negInfinity)
     965        rowLower[i] = -COIN_DBL_MAX;
     966      rowUpper[i] = LUrhs[2 * i + 1];
     967      if (rowUpper[i] >= Infinity)
     968        rowUpper[i] = COIN_DBL_MAX;
     969    }
     970    numberRows = n_con;
     971    numberColumns = n_var;
     972    numberElements = nzc;
     973    numberBinary = nbv;
     974    numberIntegers = niv;
     975    /* put in primalSolution if exists */
     976    if (X0) {
     977      primalSolution = (double *)malloc(n_var * sizeof(double));
     978      memcpy(primalSolution, X0, n_var * sizeof(double));
     979    }
     980    //double * dualSolution=NULL;
     981    if (niv + nbv > 0)
     982      mip_stuff(); // get any extra info
     983    if ((!(niv + nbv) && (csd->kind & ASL_Sufkind_input))
     984      || (rsd->kind & ASL_Sufkind_input)) {
     985      /* convert status - need info on map */
     986      static int map[] = { 1, 3, 1, 1, 2, 1, 1 };
     987      stat_map(columnStatus, n_var, map, 6, "incoming columnStatus");
     988      stat_map(rowStatus, n_con, map, 6, "incoming rowStatus");
     989    } else {
     990      /* all slack basis */
     991      // leave status for output */
     992#ifdef JJF_ZERO
     993      free(rowStatus);
     994      rowStatus = NULL;
     995      free(columnStatus);
     996      columnStatus = NULL;
     997#endif
     998    }
     999#if COIN_BIG_INDEX == 0
     1000    CoinPackedMatrix columnCopy(true, numberRows, numberColumns, numberElements,
     1001      A_vals, A_rownos, A_colstarts, NULL);
     1002#else
     1003    CoinPackedMatrix columnCopy(true, numberRows, numberColumns, numberElements,
     1004      A_vals, A_rownos,
     1005      reinterpret_cast< const CoinBigIndex * >(A_colstartsZ), NULL);
     1006#endif
     1007    matrixByRow.reverseOrderedCopyOf(columnCopy);
     1008  } else if (nonLinear == 1) {
     1009    // quadratic
     1010    asl = ASL_alloc(ASL_read_fg);
     1011    nl = jac0dim(stub, (ftnlen)strlen(stub));
     1012    free(stub);
     1013    suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
     1014    /* read  model*/
     1015    X0 = (double *)malloc(n_var * sizeof(double));
     1016    CoinZeroN(X0, n_var);
     1017    qp_read(nl, 0);
     1018    assert(n_obj == 1);
     1019    int nz = 1 + n_con;
     1020    colqp = (fint **)malloc(nz * (2 * sizeof(int *) + sizeof(double *)));
     1021    fint **rowqp = colqp + nz;
     1022    double **delsqp = (double **)(rowqp + nz);
     1023    z = (int *)malloc(nz * sizeof(int));
     1024    for (i = 0; i <= n_con; i++) {
     1025      z[i] = nqpcheck(-i, rowqp + i, colqp + i, delsqp + i);
     1026    }
     1027    qp_opify();
     1028    /* objective*/
     1029    objective = (double *)malloc(n_var * sizeof(double));
     1030    for (i = 0; i < n_var; i++)
     1031      objective[i] = 0.0;
     1032    if (n_obj) {
     1033      for (og = Ograd[0]; og; og = og->next)
     1034        objective[og->varno] = og->coef;
     1035    }
     1036    if (objtype[0])
     1037      direction = -1.0;
     1038    else
     1039      direction = 1.0;
     1040    objectiveOffset_ = objconst(0);
     1041    /* Column bounds*/
     1042    columnLower = (double *)malloc(n_var * sizeof(double));
     1043    columnUpper = (double *)malloc(n_var * sizeof(double));
     1044    for (i = 0; i < n_var; i++) {
     1045      columnLower[i] = LUv[2 * i];
     1046      if (columnLower[i] <= negInfinity)
     1047        columnLower[i] = -COIN_DBL_MAX;
     1048      columnUpper[i] = LUv[2 * i + 1];
     1049      if (columnUpper[i] >= Infinity)
     1050        columnUpper[i] = COIN_DBL_MAX;
     1051    }
     1052    // Build by row from scratch
     1053    //matrixByRow.reserve(n_var,nzc,true);
     1054    // say row orderded
     1055    matrixByRow.transpose();
     1056    /* Row bounds*/
     1057    rowLower = (double *)malloc(n_con * sizeof(double));
     1058    rowUpper = (double *)malloc(n_con * sizeof(double));
     1059    CoinBigIndex *rowStart = new CoinBigIndex[n_con + 1];
     1060    int *column = new int[nzc];
     1061    double *element = new double[nzc];
     1062    rowStart[0] = 0;
     1063    numberElements = 0;
     1064    for (i = 0; i < n_con; i++) {
     1065      rowLower[i] = LUrhs[2 * i];
     1066      if (rowLower[i] <= negInfinity)
     1067        rowLower[i] = -COIN_DBL_MAX;
     1068      rowUpper[i] = LUrhs[2 * i + 1];
     1069      if (rowUpper[i] >= Infinity)
     1070        rowUpper[i] = COIN_DBL_MAX;
     1071      for (cgrad *cg = Cgrad[i]; cg; cg = cg->next) {
     1072        column[numberElements] = cg->varno;
     1073        element[numberElements++] = cg->coef;
     1074      }
     1075      rowStart[i + 1] = numberElements;
     1076    }
     1077    assert(numberElements == nzc);
     1078    matrixByRow.appendRows(n_con, rowStart, column, element);
     1079    delete[] rowStart;
     1080    delete[] column;
     1081    delete[] element;
     1082    numberRows = n_con;
     1083    numberColumns = n_var;
     1084    //numberElements=nzc;
     1085    numberBinary = nbv;
     1086    numberIntegers = niv;
     1087    numberAllNonLinearBoth = nlvb;
     1088    numberIntegerNonLinearBoth = nlvbi;
     1089    numberAllNonLinearConstraints = nlvc;
     1090    numberIntegerNonLinearConstraints = nlvci;
     1091    numberAllNonLinearObjective = nlvo;
     1092    numberIntegerNonLinearObjective = nlvoi;
     1093    /* say we want primal solution */
     1094    want_xpi0 = 1;
     1095    //double * dualSolution=NULL;
     1096    // save asl
     1097    // Fix memory leak one day
     1098    ClpAmplInfo *info = new ClpAmplInfo;
     1099    //amplGamsData_ = info;
     1100    info->asl_ = NULL; // as wrong form asl;
     1101    info->nz_h_full_ = -1; // number of nonzeros in hessian
     1102    info->objval_called_with_current_x_ = false;
     1103    info->nerror_ = 0;
     1104    info->obj_sign_ = direction;
     1105    info->conval_called_with_current_x_ = false;
     1106    info->non_const_x_ = NULL;
     1107    info->jacval_called_with_current_x_ = false;
     1108    info->rowStart_ = NULL;
     1109    info->column_ = NULL;
     1110    info->gradient_ = NULL;
     1111    info->constraintValues_ = NULL;
     1112  } else if (nonLinear == 2) {
     1113    // General nonlinear!
     1114    //ASL_pfgh* asl = (ASL_pfgh*)ASL_alloc(ASL_read_pfgh);
     1115    asl = ASL_alloc(ASL_read_pfgh);
     1116    nl = jac0dim(stub, (ftnlen)strlen(stub));
     1117    free(stub);
     1118    suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
     1119    /* read  model*/
     1120    X0 = (double *)malloc(n_var * sizeof(double));
     1121    CoinZeroN(X0, n_var);
     1122    // code stolen from Ipopt
     1123    int retcode = pfgh_read(nl, ASL_return_read_err | ASL_findgroups);
     1124
     1125    switch (retcode) {
     1126    case ASL_readerr_none: {
     1127    } break;
     1128    case ASL_readerr_nofile: {
     1129      printf("Cannot open .nl file\n");
     1130      exit(-1);
     1131    } break;
     1132    case ASL_readerr_nonlin: {
     1133      assert(false); // this better not be an error!
     1134      printf("model involves nonlinearities (ed0read)\n");
     1135      exit(-1);
     1136    } break;
     1137    case ASL_readerr_argerr: {
     1138      printf("user-defined function with bad args\n");
     1139      exit(-1);
     1140    } break;
     1141    case ASL_readerr_unavail: {
     1142      printf("user-defined function not available\n");
     1143      exit(-1);
     1144    } break;
     1145    case ASL_readerr_corrupt: {
     1146      printf("corrupt .nl file\n");
     1147      exit(-1);
     1148    } break;
     1149    case ASL_readerr_bug: {
     1150      printf("bug in .nl reader\n");
     1151      exit(-1);
     1152    } break;
     1153    case ASL_readerr_CLP: {
     1154      printf("ASL error message: \"solver cannot handle CLP extensions\"\n");
     1155      exit(-1);
     1156    } break;
     1157    default: {
     1158      printf("Unknown error in stub file read. retcode = %d\n", retcode);
     1159      exit(-1);
     1160    } break;
     1161    }
     1162
     1163    // see "changes" in solvers directory of ampl code...
     1164    hesset(1, 0, 1, 0, nlc);
     1165
     1166    assert(n_obj == 1);
     1167    // find the nonzero structure for the hessian
     1168    // parameters to sphsetup:
     1169    int coeff_obj = 1; // coefficient of the objective fn ???
     1170    int mult_supplied = 1; // multipliers will be supplied
     1171    int uptri = 1; // only need the upper triangular part
     1172    // save asl
     1173    // Fix memory leak one day
     1174    ClpAmplInfo *info = new ClpAmplInfo;
     1175    moreInfo_ = (void *)info;
     1176    //amplGamsData_ = info;
     1177    info->asl_ = (ASL_pfgh *)asl;
     1178    // This is not easy to get from ampl so save
     1179    info->nz_h_full_ = sphsetup(-1, coeff_obj, mult_supplied, uptri);
     1180    info->objval_called_with_current_x_ = false;
     1181    info->nerror_ = 0;
     1182    info->obj_sign_ = direction;
     1183    info->conval_called_with_current_x_ = false;
     1184    info->non_const_x_ = NULL;
     1185    info->jacval_called_with_current_x_ = false;
     1186    // Look at nonlinear
     1187    if (nzc) {
     1188      n_conjac[1] = nlc; // just nonlinear
     1189      int *rowStart = new int[nlc + 1];
     1190      info->rowStart_ = rowStart;
     1191      // See how many
     1192      int current_nz = 0;
     1193      for (int i = 0; i < nlc; i++) {
     1194        for (cgrad *cg = Cgrad[i]; cg; cg = cg->next) {
     1195          current_nz++;
     1196        }
     1197      }
     1198      // setup the structure
     1199      int *column = new int[current_nz];
     1200      info->column_ = column;
     1201      current_nz = 0;
     1202      rowStart[0] = 0;
     1203      for (int i = 0; i < nlc; i++) {
     1204        for (cgrad *cg = Cgrad[i]; cg; cg = cg->next) {
     1205          cg->goff = current_nz;
     1206          //iRow[cg->goff] = i ;
     1207          //jCol[cg->goff] = cg->varno + 1;
     1208          column[cg->goff] = cg->varno;
     1209          current_nz++;
     1210        }
     1211        rowStart[i + 1] = current_nz;
     1212      }
     1213      info->gradient_ = new double[nzc];
     1214      info->constraintValues_ = new double[nlc];
     1215    }
     1216    /* objective*/
     1217    objective = (double *)malloc(n_var * sizeof(double));
     1218    for (i = 0; i < n_var; i++)
     1219      objective[i] = 0.0;
     1220    if (n_obj) {
     1221      for (og = Ograd[0]; og; og = og->next)
     1222        objective[og->varno] = og->coef;
     1223    }
     1224    if (objtype[0])
     1225      direction = -1.0;
     1226    else
     1227      direction = 1.0;
     1228    objectiveOffset_ = objconst(0);
     1229    /* Column bounds*/
     1230    columnLower = (double *)malloc(n_var * sizeof(double));
     1231    columnUpper = (double *)malloc(n_var * sizeof(double));
     1232    for (i = 0; i < n_var; i++) {
     1233      columnLower[i] = LUv[2 * i];
     1234      if (columnLower[i] <= negInfinity)
     1235        columnLower[i] = -COIN_DBL_MAX;
     1236      columnUpper[i] = LUv[2 * i + 1];
     1237      if (columnUpper[i] >= Infinity)
     1238        columnUpper[i] = COIN_DBL_MAX;
     1239    }
     1240    // Build by row from scratch
     1241    //matrixByRow.reserve(n_var,nzc,true);
     1242    // say row orderded
     1243    matrixByRow.transpose();
     1244    CoinBigIndex *rowStart = new CoinBigIndex[n_con + 1];
     1245    int *column = new int[nzc];
     1246    double *element = new double[nzc];
     1247    rowStart[0] = 0;
     1248    numberElements = 0;
     1249    /* Row bounds*/
     1250    rowLower = (double *)malloc(n_con * sizeof(double));
     1251    rowUpper = (double *)malloc(n_con * sizeof(double));
     1252    for (i = 0; i < n_con; i++) {
     1253      rowLower[i] = LUrhs[2 * i];
     1254      if (rowLower[i] <= negInfinity)
     1255        rowLower[i] = -COIN_DBL_MAX;
     1256      rowUpper[i] = LUrhs[2 * i + 1];
     1257      if (rowUpper[i] >= Infinity)
     1258        rowUpper[i] = COIN_DBL_MAX;
     1259      for (cgrad *cg = Cgrad[i]; cg; cg = cg->next) {
     1260        column[numberElements] = cg->varno;
     1261        double value = cg->coef;
     1262        if (!value)
     1263          value = -1.2345e-29;
     1264        element[numberElements++] = value;
     1265      }
     1266      rowStart[i + 1] = numberElements;
     1267    }
     1268    assert(numberElements == nzc);
     1269    matrixByRow.appendRows(n_con, rowStart, column, element);
     1270    delete[] rowStart;
     1271    delete[] column;
     1272    delete[] element;
     1273    numberRows = n_con;
     1274    numberColumns = n_var;
     1275    numberElements = nzc;
     1276    numberBinary = nbv;
     1277    numberIntegers = niv;
     1278    numberAllNonLinearBoth = nlvb;
     1279    numberIntegerNonLinearBoth = nlvbi;
     1280    numberAllNonLinearConstraints = nlvc;
     1281    numberIntegerNonLinearConstraints = nlvci;
     1282    numberAllNonLinearObjective = nlvo;
     1283    numberIntegerNonLinearObjective = nlvoi;
     1284    /* say we want primal solution */
     1285    want_xpi0 = 1;
     1286    //double * dualSolution=NULL;
     1287  } else {
     1288    abort();
     1289  }
     1290  // set problem name
     1291  problemName_ = "???";
     1292
     1293  // Build by row from scratch
     1294  const double *element = matrixByRow.getElements();
     1295  const int *column = matrixByRow.getIndices();
     1296  const CoinBigIndex *rowStart = matrixByRow.getVectorStarts();
     1297  const int *rowLength = matrixByRow.getVectorLengths();
     1298  for (i = 0; i < numberRows; i++) {
     1299    addRow(rowLength[i], column + rowStart[i],
     1300      element + rowStart[i], rowLower[i], rowUpper[i]);
     1301  }
     1302  // Now do column part
     1303  for (i = 0; i < numberColumns; i++) {
     1304    setColumnBounds(i, columnLower[i], columnUpper[i]);
     1305    setColumnObjective(i, objective[i]);
     1306  }
     1307  for (i = numberColumns - numberBinary - numberIntegers;
     1308       i < numberColumns; i++) {
     1309    setColumnIsInteger(i, true);
     1310  }
     1311  // and non linear
     1312  for (i = numberAllNonLinearBoth - numberIntegerNonLinearBoth;
     1313       i < numberAllNonLinearBoth; i++) {
     1314    setColumnIsInteger(i, true);
     1315  }
     1316  for (i = numberAllNonLinearConstraints - numberIntegerNonLinearConstraints;
     1317       i < numberAllNonLinearConstraints; i++) {
     1318    setColumnIsInteger(i, true);
     1319  }
     1320  for (i = numberAllNonLinearObjective - numberIntegerNonLinearObjective;
     1321       i < numberAllNonLinearObjective; i++) {
     1322    setColumnIsInteger(i, true);
     1323  }
     1324  free(columnLower);
     1325  free(columnUpper);
     1326  free(rowLower);
     1327  free(rowUpper);
     1328  free(objective);
     1329  // do names
     1330  int iRow;
     1331  for (iRow = 0; iRow < numberRows_; iRow++) {
     1332    char name[9];
     1333    sprintf(name, "r%7.7d", iRow);
     1334    setRowName(iRow, name);
     1335  }
     1336  int iColumn;
     1337  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     1338    char name[9];
     1339    sprintf(name, "c%7.7d", iColumn);
     1340    setColumnName(iColumn, name);
     1341  }
     1342  if (colqp) {
     1343    // add in quadratic
     1344    int nz = 1 + n_con;
     1345    int nOdd = 0;
     1346    fint **rowqp = colqp + nz;
     1347    double **delsqp = (double **)(rowqp + nz);
     1348    for (i = 0; i <= n_con; i++) {
     1349      int nels = z[i];
     1350      if (nels) {
     1351        double *element = delsqp[i];
     1352        int *start = (int *)colqp[i];
     1353        int *row = (int *)rowqp[i];
     1354        if (!element) {
     1355          // odd row - probably not quadratic
     1356          nOdd++;
     1357          continue;
     1358        }
     1359#ifdef JJF_ZERO
     1360        printf("%d quadratic els\n", nels);
     1361        for (int j = 0; j < n_var; j++) {
     1362          for (int k = start[j]; k < start[j + 1]; k++)
     1363            printf("%d %d %g\n", j, row[k], element[k]);
     1364        }
     1365#endif
     1366        if (i) {
     1367          int iRow = i - 1;
     1368          for (int j = 0; j < n_var; j++) {
     1369            for (int k = start[j]; k < start[j + 1]; k++) {
     1370              int kColumn = row[k];
     1371              double value = element[k];
     1372              // ampl gives twice with assumed 0.5
     1373              if (kColumn < j)
     1374                continue;
     1375              else if (kColumn == j)
     1376                value *= 0.5;
     1377              const char *expr = getElementAsString(iRow, j);
     1378              double constant = 0.0;
     1379              bool linear;
     1380              if (expr && strcmp(expr, "Numeric")) {
     1381                linear = false;
     1382              } else {
     1383                constant = getElement(iRow, j);
     1384                linear = true;
     1385              }
     1386              char temp[1000];
     1387              char temp2[30];
     1388              if (value == 1.0)
     1389                sprintf(temp2, "c%7.7d", kColumn);
     1390              else
     1391                sprintf(temp2, "%g*c%7.7d", value, kColumn);
     1392              if (linear) {
     1393                if (!constant)
     1394                  strcpy(temp, temp2);
     1395                else if (value > 0.0)
     1396                  sprintf(temp, "%g+%s", constant, temp2);
     1397                else
     1398                  sprintf(temp, "%g%s", constant, temp2);
     1399              } else {
     1400                if (value > 0.0)
     1401                  sprintf(temp, "%s+%s", expr, temp2);
     1402                else
     1403                  sprintf(temp, "%s%s", expr, temp2);
     1404              }
     1405              assert(strlen(temp) < 1000);
     1406              setElement(iRow, j, temp);
     1407              if (amplInfo->logLevel > 1)
     1408                printf("el for row %d column c%7.7d is %s\n", iRow, j, temp);
    4561409            }
     1410          }
     1411        } else {
     1412          // objective
     1413          for (int j = 0; j < n_var; j++) {
     1414            for (int k = start[j]; k < start[j + 1]; k++) {
     1415              int kColumn = row[k];
     1416              double value = element[k];
     1417              // ampl gives twice with assumed 0.5
     1418              if (kColumn < j)
     1419                continue;
     1420              else if (kColumn == j)
     1421                value *= 0.5;
     1422              const char *expr = getColumnObjectiveAsString(j);
     1423              double constant = 0.0;
     1424              bool linear;
     1425              if (expr && strcmp(expr, "Numeric")) {
     1426                linear = false;
     1427              } else {
     1428                constant = getColumnObjective(j);
     1429                linear = true;
     1430              }
     1431              char temp[1000];
     1432              char temp2[30];
     1433              if (value == 1.0)
     1434                sprintf(temp2, "c%7.7d", kColumn);
     1435              else
     1436                sprintf(temp2, "%g*c%7.7d", value, kColumn);
     1437              if (linear) {
     1438                if (!constant)
     1439                  strcpy(temp, temp2);
     1440                else if (value > 0.0)
     1441                  sprintf(temp, "%g+%s", constant, temp2);
     1442                else
     1443                  sprintf(temp, "%g%s", constant, temp2);
     1444              } else {
     1445                if (value > 0.0)
     1446                  sprintf(temp, "%s+%s", expr, temp2);
     1447                else
     1448                  sprintf(temp, "%s%s", expr, temp2);
     1449              }
     1450              assert(strlen(temp) < 1000);
     1451              setObjective(j, temp);
     1452              if (amplInfo->logLevel > 1)
     1453                printf("el for objective column c%7.7d is %s\n", j, temp);
     1454            }
     1455          }
    4571456        }
    458 
    459         /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
    460         Oinfo.uinfo = tempBuffer;
    461         if (getopts(argv, &Oinfo))
    462             return 1;
    463         /* objective*/
    464         obj = (double *) malloc(n_var * sizeof(double));
    465         for (i = 0; i < n_var; i++)
    466             obj[i] = 0.0;
    467         if (n_obj) {
    468             for (og = Ograd[0]; og; og = og->next)
    469                 obj[og->varno] = og->coef;
    470         }
    471         if (objtype[0])
    472             info->direction = -1.0;
    473         else
    474             info->direction = 1.0;
    475         info->offset = objconst(0);
    476         /* Column bounds*/
    477         columnLower = (double *) malloc(n_var * sizeof(double));
    478         columnUpper = (double *) malloc(n_var * sizeof(double));
    479         for (i = 0; i < n_var; i++) {
    480             columnLower[i] = LUv[2*i];
    481             if (columnLower[i] <= negInfinity)
    482                 columnLower[i] = -COIN_DBL_MAX;
    483             columnUpper[i] = LUv[2*i+1];
    484             if (columnUpper[i] >= Infinity)
    485                 columnUpper[i] = COIN_DBL_MAX;
    486         }
    487         /* Row bounds*/
    488         rowLower = (double *) malloc(n_con * sizeof(double));
    489         rowUpper = (double *) malloc(n_con * sizeof(double));
    490         for (i = 0; i < n_con; i++) {
    491             rowLower[i] = LUrhs[2*i];
    492             if (rowLower[i] <= negInfinity)
    493                 rowLower[i] = -COIN_DBL_MAX;
    494             rowUpper[i] = LUrhs[2*i+1];
    495             if (rowUpper[i] >= Infinity)
    496                 rowUpper[i] = COIN_DBL_MAX;
    497         }
    498         info->numberRows = n_con;
    499         info->numberColumns = n_var;
    500         info->numberElements = nzc;
    501         info->numberBinary = nbv;
    502         info->numberIntegers = niv + nbv;
    503         info->objective = obj;
    504         info->rowLower = rowLower;
    505         info->rowUpper = rowUpper;
    506         info->columnLower = columnLower;
    507         info->columnUpper = columnUpper;
    508 #if COIN_BIG_INDEX==0
    509         info->starts = A_colstarts;
    510 #else
    511         info->starts = A_colstartsZ;
    512 #endif
    513         /*A_colstarts=NULL;*/
    514         info->rows = A_rownos;
    515         /*A_rownos=NULL;*/
    516         info->elements = A_vals;
    517         /*A_vals=NULL;*/
    518         info->primalSolution = NULL;
    519         /* put in primalSolution if exists */
    520         if (X0) {
    521             info->primalSolution = (double *) malloc(n_var * sizeof(double));
    522             memcpy(info->primalSolution, X0, n_var*sizeof(double));
    523         }
    524         info->dualSolution = NULL;
    525         if (niv + nbv > 0)
    526             mip_stuff(); // get any extra info
    527         if ((!(niv + nbv) && (csd->kind & ASL_Sufkind_input))
    528                 || (rsd->kind & ASL_Sufkind_input)) {
    529             /* convert status - need info on map */
    530             static int map[] = {1, 3, 1, 1, 2, 1, 1};
    531             stat_map(info->columnStatus, n_var, map, 6, "incoming columnStatus");
    532             stat_map(info->rowStatus, n_con, map, 6, "incoming rowStatus");
    533         } else {
    534             /* all slack basis */
    535             // leave status for output */
    536 #ifdef JJF_ZERO
    537             free(info->rowStatus);
    538             info->rowStatus = NULL;
    539             free(info->columnStatus);
    540             info->columnStatus = NULL;
    541 #endif
    542         }
    543     } else {
    544         // QP
    545         // Add .nl if not there
    546         if (!strstr(fileName, ".nl"))
    547             strcat(fileName, ".nl");
    548         CoinModel * model = new CoinModel((nonLinearType > 10) ? 2 : 1, fileName, info);
    549         if (model->numberRows() > 0 || model->numberColumns() > 0)
    550             *coinModel = (void *) model;
    551         Oinfo.uinfo = tempBuffer;
    552         if (getopts(argv, &Oinfo))
    553             return 1;
    554         Oinfo.wantsol = 1;
    555         if (objtype[0])
    556             info->direction = -1.0;
    557         else
    558             info->direction = 1.0;
    559         model->setOptimizationDirection(info->direction);
    560         info->offset = objconst(0);
    561         info->numberRows = n_con;
    562         info->numberColumns = n_var;
    563         info->numberElements = nzc;
    564         info->numberBinary = nbv;
    565         int numberIntegers = niv + nlvci + nlvoi + nbv;
    566         if (nlvci + nlvoi + nlvc + nlvo) {
    567             // Non linear
    568             // No idea if there are overlaps so compute
    569             int numberIntegers = 0;
    570             for ( i = 0; i < n_var; i++) {
    571                 if (model->columnIsInteger(i))
    572                     numberIntegers++;
    573             }
    574         }
    575         info->numberIntegers = numberIntegers;
    576         // Say nonlinear if it is
    577         info->nonLinear = nlvc + nlvo;
    578         if (numberIntegers > 0) {
    579             mip_stuff(); // get any extra info
    580             if (info->cut)
    581                 model->setCutMarker(info->numberRows, info->cut);
    582             if (info->priorities)
    583                 model->setPriorities(info->numberColumns, info->priorities);
    584         }
    585     }
    586     /* add -solve - unless something there already
    587      - also check for sleep=yes */
    588     {
    589         int found = 0;
    590         int foundLog = 0;
    591         int foundSleep = 0;
    592         const char * something[] = {"solve", "branch", "duals", "primals", "user"};
    593         for (i = 0; i < info->numberArguments; i++) {
    594             unsigned int j;
    595             const char * argument = info->arguments[i];
    596             for (j = 0; j < sizeof(something) / sizeof(char *); j++) {
    597                 const char * check = something[j];
    598                 if (!strncmp(argument, check, sizeof(check))) {
    599                     found = (int)(j + 1);
    600                 } else if (!strncmp(argument, "log", 3)) {
    601                     foundLog = 1;
    602                 } else if (!strncmp(argument, "sleep", 5)) {
    603                     foundSleep = 1;
    604                 }
    605             }
    606         }
    607         if (foundLog) {
    608             /* print options etc */
    609             for (i = 0; i < saveArgc; i++)
    610                 printf("%s ", saveArgv[i]);
    611             printf("\n");
    612             if (environment)
    613                 printf("env %s\n", environment);
    614             /*printf("%d rows %d columns %d elements\n",n_con,n_var,nzc);*/
    615         }
    616         if (!found) {
    617             if (!strlen(algFound)) {
    618                 info->arguments = (char **) realloc(info->arguments, (info->numberArguments + 1) * sizeof(char *));
    619                 info->arguments[info->numberArguments++] = strdup("-solve");
    620             } else {
    621                 // use algorithm from keyword
    622                 info->arguments = (char **) realloc(info->arguments, (info->numberArguments + 1) * sizeof(char *));
    623                 info->arguments[info->numberArguments++] = strdup(algFound);
    624             }
    625         }
    626         if (foundSleep) {
    627             /* let user copy .nl file */
    628             fprintf(stderr, "You can copy .nl file %s for debug purposes or attach debugger\n", saveArgv[1]);
    629             fprintf(stderr, "Type q to quit, anything else to continue\n");
    630             int getChar = getc(stdin);
    631             if (getChar == 'q' || getChar == 'Q')
    632                 exit(1);
    633         }
    634     }
    635     /* add -quit */
    636     info->arguments = (char **) realloc(info->arguments, (info->numberArguments + 1) * sizeof(char *));
    637     info->arguments[info->numberArguments++] = strdup("-quit");
    638     return 0;
    639 }
    640 void freeArrays1(ampl_info * info)
    641 {
    642     free(info->objective);
    643     info->objective = NULL;
    644     free(info->rowLower);
    645     info->rowLower = NULL;
    646     free(info->rowUpper);
    647     info->rowUpper = NULL;
    648     free(info->columnLower);
    649     info->columnLower = NULL;
    650     free(info->columnUpper);
    651     info->columnUpper = NULL;
    652     /* this one not freed by ASL_free */
    653     free(info->elements);
    654     info->elements = NULL;
    655     free(info->primalSolution);
    656     info->primalSolution = NULL;
    657     free(info->dualSolution);
    658     info->dualSolution = NULL;
    659     /*free(info->rowStatus);
    660     info->rowStatus=NULL;
    661     free(info->columnStatus);
    662     info->columnStatus=NULL;*/
    663 }
    664 void freeArrays2(ampl_info * info)
    665 {
    666     free(info->primalSolution);
    667     info->primalSolution = NULL;
    668     free(info->dualSolution);
    669     info->dualSolution = NULL;
    670     free(info->rowStatus);
    671     info->rowStatus = NULL;
    672     free(info->columnStatus);
    673     info->columnStatus = NULL;
    674     free(info->priorities);
    675     info->priorities = NULL;
    676     free(info->branchDirection);
    677     info->branchDirection = NULL;
    678     free(info->pseudoDown);
    679     info->pseudoDown = NULL;
    680     free(info->pseudoUp);
    681     info->pseudoUp = NULL;
    682     free(info->sosType);
    683     info->sosType = NULL;
    684     free(info->sosPriority);
    685     info->sosPriority = NULL;
    686     free(info->sosStart);
    687     info->sosStart = NULL;
    688     free(info->sosIndices);
    689     info->sosIndices = NULL;
    690     free(info->sosReference);
    691     info->sosReference = NULL;
    692     free(info->cut);
    693     info->cut = NULL;
    694     ASL_free(&asl);
    695 }
    696 void freeArgs(ampl_info * info)
    697 {
    698     int i;
    699     for ( i = 0; i < info->numberArguments; i++)
    700         free(info->arguments[i]);
    701     free(info->arguments);
    702 }
    703 int ampl_obj_prec()
    704 {
    705     int precision = obj_prec();
    706     if (precision<=0)
    707         precision=15;
    708     return precision;
    709 }
    710 void writeAmpl(ampl_info * info)
    711 {
    712     char buf[1000];
    713     typedef struct {
    714         const char *msg;
    715         int code;
    716         int wantObj;
    717     } Sol_info;
    718     static Sol_info solinfo[] = {
    719         { "optimal solution",                   000, 1 },
    720         { "infeasible",                         200, 1 },
    721         { "unbounded",                          300, 0 },
    722         { "iteration limit etc",                        400, 1 },
    723         { "solution limit",                             401, 1 },
    724         { "ran out of space",                   500, 0 },
    725         { "status unknown",                             501, 1 },
    726         { "bug!",                                       502, 0 },
    727         { "best MIP solution so far restored",  101, 1 },
    728         { "failed to restore best MIP solution",        503, 1 },
    729         { "optimal (?) solution",                       100, 1 }
    730     };
    731     /* convert status - need info on map */
    732     static int map[] = {0, 3, 4, 1};
    733     sprintf(buf, "%s %s", Oinfo.bsname, info->buffer);
    734     solve_result_num = solinfo[info->problemStatus].code;
    735     if (info->columnStatus) {
    736         stat_map(info->columnStatus, n_var, map, 4, "outgoing columnStatus");
    737         stat_map(info->rowStatus, n_con, map, 4, "outgoing rowStatus");
    738         suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
    739         suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
    740     }
    741     write_sol(buf, info->primalSolution, info->dualSolution, &Oinfo);
    742 }
    743 /* Read a problem from AMPL nl file
    744  */
    745 CoinModel::CoinModel( int nonLinear, const char * fileName, const void * info)
    746         :  CoinBaseModel(),
    747         maximumRows_(0),
    748         maximumColumns_(0),
    749         numberElements_(0),
    750         maximumElements_(0),
    751         numberQuadraticElements_(0),
    752         maximumQuadraticElements_(0),
    753         rowLower_(NULL),
    754         rowUpper_(NULL),
    755         rowType_(NULL),
    756         objective_(NULL),
    757         columnLower_(NULL),
    758         columnUpper_(NULL),
    759         integerType_(NULL),
    760         columnType_(NULL),
    761         start_(NULL),
    762         elements_(NULL),
    763         packedMatrix_(NULL),
    764         quadraticElements_(NULL),
    765         sortIndices_(NULL),
    766         sortElements_(NULL),
    767         sortSize_(0),
    768         sizeAssociated_(0),
    769         associated_(NULL),
    770         numberSOS_(0),
    771         startSOS_(NULL),
    772         memberSOS_(NULL),
    773         typeSOS_(NULL),
    774         prioritySOS_(NULL),
    775         referenceSOS_(NULL),
    776         priority_(NULL),
    777         cut_(NULL),
    778         moreInfo_(NULL),
    779         type_(-1),
    780         noNames_(false),
    781         links_(0)
    782 {
    783     problemName_ = "";
    784     int status = 0;
    785     if (!strcmp(fileName, "-") || !strcmp(fileName, "stdin")) {
    786         // stdin
    787     } else {
    788         std::string name = fileName;
    789         bool readable = fileCoinReadable(name);
    790         if (!readable) {
    791             std::cerr << "Unable to open file "
    792                       << fileName << std::endl;
    793             status = -1;
    794         }
    795     }
    796     if (!status) {
    797         gdb(nonLinear, fileName, info);
    798     }
    799 }
    800 #ifdef JJF_ZERO
    801 static real
    802 qterm(ASL *asl, fint *colq, fint *rowq, real *delsq)
    803 {
    804     double t, t1, *x, *x0, *xe;
    805     fint *rq0, *rqe;
    806 
    807     t = 0.;
    808     x0 = x = X0;
    809     xe = x + n_var;
    810     rq0 = rowq;
    811     while (x < xe) {
    812         t1 = *x++;
    813         rqe = rq0 + *++colq;
    814         while (rowq < rqe)
    815             t += t1 * x0[*rowq++]**delsq++;
    816     }
    817     return 0.5 * t;
    818 }
    819 #endif
    820 // stolen from IPopt with changes
    821 typedef struct {
    822     double obj_sign_;
    823     ASL_pfgh * asl_;
    824     double * non_const_x_;
    825     int * column_; // for jacobian
    826     int * rowStart_;
    827     double * gradient_;
    828     double * constraintValues_;
    829     int nz_h_full_; // number of nonzeros in hessian
    830     int nerror_;
    831     bool objval_called_with_current_x_;
    832     bool conval_called_with_current_x_;
    833     bool jacval_called_with_current_x_;
    834 } ClpAmplInfo;
    835 
    836 void
    837 CoinModel::gdb( int nonLinear, const char * fileName, const void * info)
    838 {
    839     const ampl_info * amplInfo = (const ampl_info *) info;
    840     ograd *og = NULL;
    841     int i;
    842     SufDesc *csd = NULL;
    843     SufDesc *rsd = NULL;
    844     /*bool *basis, *lower;*/
    845     /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
    846     //char tempBuffer[20];
    847     double * objective = NULL;
    848     double * columnLower = NULL;
    849     double * columnUpper = NULL;
    850     double * rowLower = NULL;
    851     double * rowUpper = NULL;
    852     int * columnStatus = NULL;
    853     int * rowStatus = NULL;
    854     int numberRows = -1;
    855     int numberColumns = -1;
    856     int numberElements = -1;
    857     int numberBinary = -1;
    858     int numberIntegers = -1;
    859     int numberAllNonLinearBoth = 0;
    860     int numberIntegerNonLinearBoth = 0;
    861     int numberAllNonLinearConstraints = 0;
    862     int numberIntegerNonLinearConstraints = 0;
    863     int numberAllNonLinearObjective = 0;
    864     int numberIntegerNonLinearObjective = 0;
    865     double * primalSolution = NULL;
    866     double direction = 1.0;
    867     char * stub = strdup(fileName);
    868     CoinPackedMatrix matrixByRow;
    869     fint ** colqp = NULL;
    870     int *z = NULL;
    871     if (nonLinear == 0) {
    872         // linear
    873         asl = ASL_alloc(ASL_read_f);
    874         nl = jac0dim(stub, 0);
    875         free(stub);
    876         suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
    877 
    878         /* set A_vals to get the constraints column-wise (malloc so can be freed) */
    879         A_vals = (double *) malloc(nzc * sizeof(double));
    880         if (!A_vals) {
    881             printf("no memory\n");
    882             return ;
    883         }
    884         /* say we want primal solution */
    885         want_xpi0 = 1;
    886         /* for basis info */
    887         columnStatus = (int *) malloc(n_var * sizeof(int));
    888         rowStatus = (int *) malloc(n_con * sizeof(int));
    889         csd = suf_iput("sstatus", ASL_Sufkind_var, columnStatus);
    890         rsd = suf_iput("sstatus", ASL_Sufkind_con, rowStatus);
    891         /* read linear model*/
    892 #if COIN_BIG_INDEX==2
    893         f_read(nl, ASL_allow_Z|ASL_use_Z);
    894 #else
    895         f_read(nl, 0);
    896 #endif
    897         // see if any sos
    898         if (true) {
    899             char *sostype;
    900             int nsosnz, *sosbeg, *sosind, * sospri;
    901             double *sosref;
    902             int nsos;
    903             int i = ASL_suf_sos_explict_free;
    904             int copri[2], **p_sospri;
    905             copri[0] = 0;
    906             copri[1] = 0;
    907             p_sospri = &sospri;
    908             nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
    909                            &sosbeg, &sosind, &sosref);
    910             if (nsos) {
    911                 abort();
    912 #ifdef JJF_ZERO
    913                 info->numberSos = nsos;
    914                 info->sosType = (char *) malloc(nsos);
    915                 info->sosPriority = (int *) malloc(nsos * sizeof(int));
    916                 info->sosStart = (int *) malloc((nsos + 1) * sizeof(int));
    917                 info->sosIndices = (int *) malloc(nsosnz * sizeof(int));
    918                 info->sosReference = (double *) malloc(nsosnz * sizeof(double));
    919                 sos_kludge(nsos, sosbeg, sosref, sosind);
    920                 for (int i = 0; i < nsos; i++) {
    921                     int ichar = sostype[i];
    922                     assert (ichar == '1' || ichar == '2');
    923                     info->sosType[i] = ichar - '0';
    924                 }
    925                 memcpy(info->sosPriority, sospri, nsos*sizeof(int));
    926                 memcpy(info->sosStart, sosbeg, (nsos + 1)*sizeof(int));
    927                 memcpy(info->sosIndices, sosind, nsosnz*sizeof(int));
    928                 memcpy(info->sosReference, sosref, nsosnz*sizeof(double));
    929 #endif
    930             }
    931         }
    932 
    933         /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
    934         //Oinfo.uinfo = tempBuffer;
    935         //if (getopts(argv, &Oinfo))
    936         //return 1;
    937         /* objective*/
    938         objective = (double *) malloc(n_var * sizeof(double));
    939         for (i = 0; i < n_var; i++)
    940             objective[i] = 0.0;
    941         if (n_obj) {
    942             for (og = Ograd[0]; og; og = og->next)
    943                 objective[og->varno] = og->coef;
    944         }
    945         if (objtype[0])
    946             direction = -1.0;
    947         else
    948             direction = 1.0;
    949         objectiveOffset_ = objconst(0);
    950         /* Column bounds*/
    951         columnLower = (double *) malloc(n_var * sizeof(double));
    952         columnUpper = (double *) malloc(n_var * sizeof(double));
    953         for (i = 0; i < n_var; i++) {
    954             columnLower[i] = LUv[2*i];
    955             if (columnLower[i] <= negInfinity)
    956                 columnLower[i] = -COIN_DBL_MAX;
    957             columnUpper[i] = LUv[2*i+1];
    958             if (columnUpper[i] >= Infinity)
    959                 columnUpper[i] = COIN_DBL_MAX;
    960         }
    961         /* Row bounds*/
    962         rowLower = (double *) malloc(n_con * sizeof(double));
    963         rowUpper = (double *) malloc(n_con * sizeof(double));
    964         for (i = 0; i < n_con; i++) {
    965             rowLower[i] = LUrhs[2*i];
    966             if (rowLower[i] <= negInfinity)
    967                 rowLower[i] = -COIN_DBL_MAX;
    968             rowUpper[i] = LUrhs[2*i+1];
    969             if (rowUpper[i] >= Infinity)
    970                 rowUpper[i] = COIN_DBL_MAX;
    971         }
    972         numberRows = n_con;
    973         numberColumns = n_var;
    974         numberElements = nzc;
    975         numberBinary = nbv;
    976         numberIntegers = niv;
    977         /* put in primalSolution if exists */
    978         if (X0) {
    979             primalSolution = (double *) malloc(n_var * sizeof(double));
    980             memcpy( primalSolution, X0, n_var*sizeof(double));
    981         }
    982         //double * dualSolution=NULL;
    983         if (niv + nbv > 0)
    984             mip_stuff(); // get any extra info
    985         if ((!(niv + nbv) && (csd->kind & ASL_Sufkind_input))
    986                 || (rsd->kind & ASL_Sufkind_input)) {
    987             /* convert status - need info on map */
    988             static int map[] = {1, 3, 1, 1, 2, 1, 1};
    989             stat_map(columnStatus, n_var, map, 6, "incoming columnStatus");
    990             stat_map(rowStatus, n_con, map, 6, "incoming rowStatus");
    991         } else {
    992             /* all slack basis */
    993             // leave status for output */
    994 #ifdef JJF_ZERO
    995             free(rowStatus);
    996             rowStatus = NULL;
    997             free(columnStatus);
    998             columnStatus = NULL;
    999 #endif
    1000         }
    1001 #if COIN_BIG_INDEX==0
    1002         CoinPackedMatrix columnCopy(true, numberRows, numberColumns, numberElements,
    1003                                     A_vals, A_rownos, A_colstarts, NULL);
    1004 #else
    1005         CoinPackedMatrix columnCopy(true, numberRows, numberColumns, numberElements,
    1006                                     A_vals, A_rownos,
    1007                                     reinterpret_cast<const CoinBigIndex *>(A_colstartsZ), NULL);
    1008 #endif
    1009         matrixByRow.reverseOrderedCopyOf(columnCopy);
    1010     } else if (nonLinear == 1) {
    1011         // quadratic
    1012         asl = ASL_alloc(ASL_read_fg);
    1013         nl = jac0dim(stub, (ftnlen) strlen(stub));
    1014         free(stub);
    1015         suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
    1016         /* read  model*/
    1017         X0 = (double*) malloc(n_var * sizeof(double));
    1018         CoinZeroN(X0, n_var);
    1019         qp_read(nl, 0);
    1020         assert (n_obj == 1);
    1021         int nz = 1 + n_con;
    1022         colqp = (fint**) malloc(nz * (2 * sizeof(int*)
    1023                                       + sizeof(double*)));
    1024         fint ** rowqp = colqp + nz;
    1025         double ** delsqp = (double **)(rowqp + nz);
    1026         z = (int*) malloc(nz * sizeof(int));
    1027         for (i = 0; i <= n_con; i++) {
    1028             z[i] = nqpcheck(-i, rowqp + i, colqp + i, delsqp + i);
    1029         }
    1030         qp_opify();
    1031         /* objective*/
    1032         objective = (double *) malloc(n_var * sizeof(double));
    1033         for (i = 0; i < n_var; i++)
    1034             objective[i] = 0.0;
    1035         if (n_obj) {
    1036             for (og = Ograd[0]; og; og = og->next)
    1037                 objective[og->varno] = og->coef;
    1038         }
    1039         if (objtype[0])
    1040             direction = -1.0;
    1041         else
    1042             direction = 1.0;
    1043         objectiveOffset_ = objconst(0);
    1044         /* Column bounds*/
    1045         columnLower = (double *) malloc(n_var * sizeof(double));
    1046         columnUpper = (double *) malloc(n_var * sizeof(double));
    1047         for (i = 0; i < n_var; i++) {
    1048             columnLower[i] = LUv[2*i];
    1049             if (columnLower[i] <= negInfinity)
    1050                 columnLower[i] = -COIN_DBL_MAX;
    1051             columnUpper[i] = LUv[2*i+1];
    1052             if (columnUpper[i] >= Infinity)
    1053                 columnUpper[i] = COIN_DBL_MAX;
    1054         }
    1055         // Build by row from scratch
    1056         //matrixByRow.reserve(n_var,nzc,true);
    1057         // say row orderded
    1058         matrixByRow.transpose();
    1059         /* Row bounds*/
    1060         rowLower = (double *) malloc(n_con * sizeof(double));
    1061         rowUpper = (double *) malloc(n_con * sizeof(double));
    1062         CoinBigIndex * rowStart = new CoinBigIndex [n_con+1];
    1063         int * column = new int [nzc];
    1064         double * element = new double [nzc];
    1065         rowStart[0] = 0;
    1066         numberElements = 0;
    1067         for (i = 0; i < n_con; i++) {
    1068             rowLower[i] = LUrhs[2*i];
    1069             if (rowLower[i] <= negInfinity)
    1070                 rowLower[i] = -COIN_DBL_MAX;
    1071             rowUpper[i] = LUrhs[2*i+1];
    1072             if (rowUpper[i] >= Infinity)
    1073                 rowUpper[i] = COIN_DBL_MAX;
    1074             for (cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
    1075                 column[numberElements] = cg->varno;
    1076                 element[numberElements++] = cg->coef;
    1077             }
    1078             rowStart[i+1] = numberElements;
    1079         }
    1080         assert (numberElements == nzc);
    1081         matrixByRow.appendRows(n_con, rowStart, column, element);
    1082         delete [] rowStart;
    1083         delete [] column;
    1084         delete [] element;
    1085         numberRows = n_con;
    1086         numberColumns = n_var;
    1087         //numberElements=nzc;
    1088         numberBinary = nbv;
    1089         numberIntegers = niv;
    1090         numberAllNonLinearBoth = nlvb;
    1091         numberIntegerNonLinearBoth = nlvbi;
    1092         numberAllNonLinearConstraints = nlvc;
    1093         numberIntegerNonLinearConstraints = nlvci;
    1094         numberAllNonLinearObjective = nlvo;
    1095         numberIntegerNonLinearObjective = nlvoi;
    1096         /* say we want primal solution */
    1097         want_xpi0 = 1;
    1098         //double * dualSolution=NULL;
    1099         // save asl
    1100         // Fix memory leak one day
    1101         ClpAmplInfo * info = new ClpAmplInfo;
    1102         //amplGamsData_ = info;
    1103         info->asl_ = NULL; // as wrong form asl;
    1104         info->nz_h_full_ = -1; // number of nonzeros in hessian
    1105         info->objval_called_with_current_x_ = false;
    1106         info->nerror_ = 0;
    1107         info->obj_sign_ = direction;
    1108         info->conval_called_with_current_x_ = false;
    1109         info->non_const_x_ = NULL;
    1110         info->jacval_called_with_current_x_ = false;
    1111         info->rowStart_ = NULL;
    1112         info->column_ = NULL;
    1113         info->gradient_ = NULL;
    1114         info->constraintValues_ = NULL;
    1115     } else if (nonLinear == 2) {
    1116         // General nonlinear!
    1117         //ASL_pfgh* asl = (ASL_pfgh*)ASL_alloc(ASL_read_pfgh);
    1118         asl = ASL_alloc(ASL_read_pfgh);
    1119         nl = jac0dim(stub, (ftnlen) strlen(stub));
    1120         free(stub);
    1121         suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
    1122         /* read  model*/
    1123         X0 = (double*) malloc(n_var * sizeof(double));
    1124         CoinZeroN(X0, n_var);
    1125         // code stolen from Ipopt
    1126         int retcode = pfgh_read(nl, ASL_return_read_err | ASL_findgroups);
    1127 
    1128         switch (retcode) {
    1129         case ASL_readerr_none : {}
    1130         break;
    1131         case ASL_readerr_nofile : {
    1132             printf( "Cannot open .nl file\n");
    1133             exit(-1);
    1134         }
    1135         break;
    1136         case ASL_readerr_nonlin : {
    1137             assert(false); // this better not be an error!
    1138             printf( "model involves nonlinearities (ed0read)\n");
    1139             exit(-1);
    1140         }
    1141         break;
    1142         case  ASL_readerr_argerr : {
    1143             printf( "user-defined function with bad args\n");
    1144             exit(-1);
    1145         }
    1146         break;
    1147         case ASL_readerr_unavail : {
    1148             printf( "user-defined function not available\n");
    1149             exit(-1);
    1150         }
    1151         break;
    1152         case ASL_readerr_corrupt : {
    1153             printf( "corrupt .nl file\n");
    1154             exit(-1);
    1155         }
    1156         break;
    1157         case ASL_readerr_bug : {
    1158             printf( "bug in .nl reader\n");
    1159             exit(-1);
    1160         }
    1161         break;
    1162         case ASL_readerr_CLP : {
    1163             printf( "ASL error message: \"solver cannot handle CLP extensions\"\n");
    1164             exit(-1);
    1165         }
    1166         break;
    1167         default: {
    1168             printf( "Unknown error in stub file read. retcode = %d\n", retcode);
    1169             exit(-1);
    1170         }
    1171         break;
    1172         }
    1173 
    1174         // see "changes" in solvers directory of ampl code...
    1175         hesset(1, 0, 1, 0, nlc);
    1176 
    1177         assert (n_obj == 1);
    1178         // find the nonzero structure for the hessian
    1179         // parameters to sphsetup:
    1180         int coeff_obj = 1; // coefficient of the objective fn ???
    1181         int mult_supplied = 1; // multipliers will be supplied
    1182         int uptri = 1; // only need the upper triangular part
    1183         // save asl
    1184         // Fix memory leak one day
    1185         ClpAmplInfo * info = new ClpAmplInfo;
    1186         moreInfo_ = (void *) info;
    1187         //amplGamsData_ = info;
    1188         info->asl_ = (ASL_pfgh *) asl;
    1189         // This is not easy to get from ampl so save
    1190         info->nz_h_full_ = sphsetup(-1, coeff_obj, mult_supplied, uptri);
    1191         info->objval_called_with_current_x_ = false;
    1192         info->nerror_ = 0;
    1193         info->obj_sign_ = direction;
    1194         info->conval_called_with_current_x_ = false;
    1195         info->non_const_x_ = NULL;
    1196         info->jacval_called_with_current_x_ = false;
    1197         // Look at nonlinear
    1198         if (nzc) {
    1199             n_conjac[1] = nlc; // just nonlinear
    1200             int * rowStart = new int [nlc+1];
    1201             info->rowStart_ = rowStart;
    1202             // See how many
    1203             int  current_nz = 0;
    1204             for (int i = 0; i < nlc; i++) {
    1205                 for (cgrad* cg = Cgrad[i]; cg; cg = cg->next) {
    1206                     current_nz++;
    1207                 }
    1208             }
    1209             // setup the structure
    1210             int * column = new int [current_nz];
    1211             info->column_ = column;
    1212             current_nz = 0;
    1213             rowStart[0] = 0;
    1214             for (int i = 0; i < nlc; i++) {
    1215                 for (cgrad* cg = Cgrad[i]; cg; cg = cg->next) {
    1216                     cg->goff = current_nz;
    1217                     //iRow[cg->goff] = i ;
    1218                     //jCol[cg->goff] = cg->varno + 1;
    1219                     column[cg->goff] = cg->varno ;
    1220                     current_nz++;
    1221                 }
    1222                 rowStart[i+1] = current_nz;
    1223             }
    1224             info->gradient_ = new double [nzc];
    1225             info->constraintValues_ = new double [nlc];
    1226         }
    1227         /* objective*/
    1228         objective = (double *) malloc(n_var * sizeof(double));
    1229         for (i = 0; i < n_var; i++)
    1230             objective[i] = 0.0;
    1231         if (n_obj) {
    1232             for (og = Ograd[0]; og; og = og->next)
    1233                 objective[og->varno] = og->coef;
    1234         }
    1235         if (objtype[0])
    1236             direction = -1.0;
    1237         else
    1238             direction = 1.0;
    1239         objectiveOffset_ = objconst(0);
    1240         /* Column bounds*/
    1241         columnLower = (double *) malloc(n_var * sizeof(double));
    1242         columnUpper = (double *) malloc(n_var * sizeof(double));
    1243         for (i = 0; i < n_var; i++) {
    1244             columnLower[i] = LUv[2*i];
    1245             if (columnLower[i] <= negInfinity)
    1246                 columnLower[i] = -COIN_DBL_MAX;
    1247             columnUpper[i] = LUv[2*i+1];
    1248             if (columnUpper[i] >= Infinity)
    1249                 columnUpper[i] = COIN_DBL_MAX;
    1250         }
    1251         // Build by row from scratch
    1252         //matrixByRow.reserve(n_var,nzc,true);
    1253         // say row orderded
    1254         matrixByRow.transpose();
    1255         CoinBigIndex * rowStart = new CoinBigIndex [n_con+1];
    1256         int * column = new int [nzc];
    1257         double * element = new double [nzc];
    1258         rowStart[0] = 0;
    1259         numberElements = 0;
    1260         /* Row bounds*/
    1261         rowLower = (double *) malloc(n_con * sizeof(double));
    1262         rowUpper = (double *) malloc(n_con * sizeof(double));
    1263         for (i = 0; i < n_con; i++) {
    1264             rowLower[i] = LUrhs[2*i];
    1265             if (rowLower[i] <= negInfinity)
    1266                 rowLower[i] = -COIN_DBL_MAX;
    1267             rowUpper[i] = LUrhs[2*i+1];
    1268             if (rowUpper[i] >= Infinity)
    1269                 rowUpper[i] = COIN_DBL_MAX;
    1270             for (cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
    1271                 column[numberElements] = cg->varno;
    1272                 double value = cg->coef;
    1273                 if (!value)
    1274                     value = -1.2345e-29;
    1275                 element[numberElements++] = value;
    1276             }
    1277             rowStart[i+1] = numberElements;
    1278         }
    1279         assert (numberElements == nzc);
    1280         matrixByRow.appendRows(n_con, rowStart, column, element);
    1281         delete [] rowStart;
    1282         delete [] column;
    1283         delete [] element;
    1284         numberRows = n_con;
    1285         numberColumns = n_var;
    1286         numberElements = nzc;
    1287         numberBinary = nbv;
    1288         numberIntegers = niv;
    1289         numberAllNonLinearBoth = nlvb;
    1290         numberIntegerNonLinearBoth = nlvbi;
    1291         numberAllNonLinearConstraints = nlvc;
    1292         numberIntegerNonLinearConstraints = nlvci;
    1293         numberAllNonLinearObjective = nlvo;
    1294         numberIntegerNonLinearObjective = nlvoi;
    1295         /* say we want primal solution */
    1296         want_xpi0 = 1;
    1297         //double * dualSolution=NULL;
    1298     } else {
    1299         abort();
    1300     }
    1301     // set problem name
    1302     problemName_ = "???";
    1303 
    1304     // Build by row from scratch
    1305     const double * element = matrixByRow.getElements();
    1306     const int * column = matrixByRow.getIndices();
    1307     const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
    1308     const int * rowLength = matrixByRow.getVectorLengths();
    1309     for (i = 0; i < numberRows; i++) {
    1310         addRow(rowLength[i], column + rowStart[i],
    1311                element + rowStart[i], rowLower[i], rowUpper[i]);
    1312     }
    1313     // Now do column part
    1314     for (i = 0; i < numberColumns; i++) {
    1315         setColumnBounds(i, columnLower[i], columnUpper[i]);
    1316         setColumnObjective(i, objective[i]);
    1317     }
    1318     for ( i = numberColumns - numberBinary - numberIntegers;
    1319             i < numberColumns; i++) {
    1320         setColumnIsInteger(i, true);
    1321     }
    1322     // and non linear
    1323     for (i = numberAllNonLinearBoth - numberIntegerNonLinearBoth;
    1324             i < numberAllNonLinearBoth; i++) {
    1325         setColumnIsInteger(i, true);
    1326     }
    1327     for (i = numberAllNonLinearConstraints - numberIntegerNonLinearConstraints;
    1328             i < numberAllNonLinearConstraints; i++) {
    1329         setColumnIsInteger(i, true);
    1330     }
    1331     for (i = numberAllNonLinearObjective - numberIntegerNonLinearObjective;
    1332             i < numberAllNonLinearObjective; i++) {
    1333         setColumnIsInteger(i, true);
    1334     }
    1335     free(columnLower);
    1336     free(columnUpper);
    1337     free(rowLower);
    1338     free(rowUpper);
    1339     free(objective);
    1340     // do names
    1341     int iRow;
    1342     for (iRow = 0; iRow < numberRows_; iRow++) {
    1343         char name[9];
    1344         sprintf(name, "r%7.7d", iRow);
    1345         setRowName(iRow, name);
    1346     }
    1347     int iColumn;
    1348     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1349         char name[9];
    1350         sprintf(name, "c%7.7d", iColumn);
    1351         setColumnName(iColumn, name);
    1352     }
    1353     if (colqp) {
    1354         // add in quadratic
    1355         int nz = 1 + n_con;
    1356         int nOdd = 0;
    1357         fint ** rowqp = colqp + nz;
    1358         double ** delsqp = (double **)(rowqp + nz);
    1359         for (i = 0; i <= n_con; i++) {
    1360             int nels = z[i];
    1361             if (nels) {
    1362                 double * element = delsqp[i];
    1363                 int * start = (int *) colqp[i];
    1364                 int * row = (int *) rowqp[i];
    1365                 if (!element) {
    1366                     // odd row - probably not quadratic
    1367                     nOdd++;
    1368                     continue;
    1369                 }
    1370 #ifdef JJF_ZERO
    1371                 printf("%d quadratic els\n", nels);
    1372                 for (int j = 0; j < n_var; j++) {
    1373                     for (int k = start[j]; k < start[j+1]; k++)
    1374                         printf("%d %d %g\n", j, row[k], element[k]);
    1375                 }
    1376 #endif
    1377                 if (i) {
    1378                     int iRow = i - 1;
    1379                     for (int j = 0; j < n_var; j++) {
    1380                         for (int k = start[j]; k < start[j+1]; k++) {
    1381                             int kColumn = row[k];
    1382                             double value = element[k];
    1383                             // ampl gives twice with assumed 0.5
    1384                             if (kColumn < j)
    1385                                 continue;
    1386                             else if (kColumn == j)
    1387                                 value *= 0.5;
    1388                             const char * expr = getElementAsString(iRow, j);
    1389                             double constant = 0.0;
    1390                             bool linear;
    1391                             if (expr && strcmp(expr, "Numeric")) {
    1392                                 linear = false;
    1393                             } else {
    1394                                 constant = getElement(iRow, j);
    1395                                 linear = true;
    1396                             }
    1397                             char temp[1000];
    1398                             char temp2[30];
    1399                             if (value == 1.0)
    1400                                 sprintf(temp2, "c%7.7d", kColumn);
    1401                             else
    1402                                 sprintf(temp2, "%g*c%7.7d", value, kColumn);
    1403                             if (linear) {
    1404                                 if (!constant)
    1405                                     strcpy(temp, temp2);
    1406                                 else if (value > 0.0)
    1407                                     sprintf(temp, "%g+%s", constant, temp2);
    1408                                 else
    1409                                     sprintf(temp, "%g%s", constant, temp2);
    1410                             } else {
    1411                                 if (value > 0.0)
    1412                                     sprintf(temp, "%s+%s", expr, temp2);
    1413                                 else
    1414                                     sprintf(temp, "%s%s", expr, temp2);
    1415                             }
    1416                             assert (strlen(temp) < 1000);
    1417                             setElement(iRow, j, temp);
    1418                             if (amplInfo->logLevel > 1)
    1419                                 printf("el for row %d column c%7.7d is %s\n", iRow, j, temp);
    1420                         }
    1421                     }
    1422                 } else {
    1423                     // objective
    1424                     for (int j = 0; j < n_var; j++) {
    1425                         for (int k = start[j]; k < start[j+1]; k++) {
    1426                             int kColumn = row[k];
    1427                             double value = element[k];
    1428                             // ampl gives twice with assumed 0.5
    1429                             if (kColumn < j)
    1430                                 continue;
    1431                             else if (kColumn == j)
    1432                                 value *= 0.5;
    1433                             const char * expr = getColumnObjectiveAsString(j);
    1434                             double constant = 0.0;
    1435                             bool linear;
    1436                             if (expr && strcmp(expr, "Numeric")) {
    1437                                 linear = false;
    1438                             } else {
    1439                                 constant = getColumnObjective(j);
    1440                                 linear = true;
    1441                             }
    1442                             char temp[1000];
    1443                             char temp2[30];
    1444                             if (value == 1.0)
    1445                                 sprintf(temp2, "c%7.7d", kColumn);
    1446                             else
    1447                                 sprintf(temp2, "%g*c%7.7d", value, kColumn);
    1448                             if (linear) {
    1449                                 if (!constant)
    1450                                     strcpy(temp, temp2);
    1451                                 else if (value > 0.0)
    1452                                     sprintf(temp, "%g+%s", constant, temp2);
    1453                                 else
    1454                                     sprintf(temp, "%g%s", constant, temp2);
    1455                             } else {
    1456                                 if (value > 0.0)
    1457                                     sprintf(temp, "%s+%s", expr, temp2);
    1458                                 else
    1459                                     sprintf(temp, "%s%s", expr, temp2);
    1460                             }
    1461                             assert (strlen(temp) < 1000);
    1462                             setObjective(j, temp);
    1463                             if (amplInfo->logLevel > 1)
    1464                                 printf("el for objective column c%7.7d is %s\n", j, temp);
    1465                         }
    1466                     }
    1467                 }
    1468             }
    1469         }
    1470         if (nOdd) {
    1471             printf("%d non-linear constraints could not be converted to quadratic\n", nOdd);
    1472             exit(77);
    1473         }
    1474     }
    1475     free(colqp);
    1476     free(z);
    1477     // see if any sos
    1478     {
    1479         char *sostype;
    1480         int nsosnz, *sosbeg, *sosind, * sospri;
    1481         double *sosref;
    1482         int nsos;
    1483         int i = ASL_suf_sos_explict_free;
    1484         int copri[2], **p_sospri;
    1485         copri[0] = 0;
    1486         copri[1] = 0;
    1487         p_sospri = &sospri;
    1488         nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
    1489                        &sosbeg, &sosind, &sosref);
    1490         if (nsos) {
    1491             numberSOS_ = nsos;
    1492             typeSOS_ = new int [numberSOS_];
    1493             prioritySOS_ = new int [numberSOS_];
    1494             startSOS_ = new int [numberSOS_+1];
    1495             memberSOS_ = new int[nsosnz];
    1496             referenceSOS_ = new double [nsosnz];
    1497             sos_kludge(nsos, sosbeg, sosref, sosind);
    1498             for (int i = 0; i < nsos; i++) {
    1499                 int ichar = sostype[i];
    1500                 assert (ichar == '1' || ichar == '2');
    1501                 typeSOS_[i] = ichar - '0';
    1502             }
    1503             memcpy(prioritySOS_, sospri, nsos*sizeof(int));
    1504             memcpy(startSOS_, sosbeg, (nsos + 1)*sizeof(int));
    1505             memcpy(memberSOS_, sosind, nsosnz*sizeof(int));
    1506             memcpy(referenceSOS_, sosref, nsosnz*sizeof(double));
    1507         }
    1508     }
     1457      }
     1458    }
     1459    if (nOdd) {
     1460      printf("%d non-linear constraints could not be converted to quadratic\n", nOdd);
     1461      exit(77);
     1462    }
     1463  }
     1464  free(colqp);
     1465  free(z);
     1466  // see if any sos
     1467  {
     1468    char *sostype;
     1469    int nsosnz, *sosbeg, *sosind, *sospri;
     1470    double *sosref;
     1471    int nsos;
     1472    int i = ASL_suf_sos_explict_free;
     1473    int copri[2], **p_sospri;
     1474    copri[0] = 0;
     1475    copri[1] = 0;
     1476    p_sospri = &sospri;
     1477    nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
     1478      &sosbeg, &sosind, &sosref);
     1479    if (nsos) {
     1480      numberSOS_ = nsos;
     1481      typeSOS_ = new int[numberSOS_];
     1482      prioritySOS_ = new int[numberSOS_];
     1483      startSOS_ = new int[numberSOS_ + 1];
     1484      memberSOS_ = new int[nsosnz];
     1485      referenceSOS_ = new double[nsosnz];
     1486      sos_kludge(nsos, sosbeg, sosref, sosind);
     1487      for (int i = 0; i < nsos; i++) {
     1488        int ichar = sostype[i];
     1489        assert(ichar == '1' || ichar == '2');
     1490        typeSOS_[i] = ichar - '0';
     1491      }
     1492      memcpy(prioritySOS_, sospri, nsos * sizeof(int));
     1493      memcpy(startSOS_, sosbeg, (nsos + 1) * sizeof(int));
     1494      memcpy(memberSOS_, sosind, nsosnz * sizeof(int));
     1495      memcpy(referenceSOS_, sosref, nsosnz * sizeof(double));
     1496    }
     1497  }
    15091498}
    15101499#else
    15111500#include "Clp_ampl.h"
    1512 int
    1513 readAmpl(ampl_info * , int , char **, void ** )
    1514 {
    1515     return 0;
     1501int readAmpl(ampl_info *, int, char **, void **)
     1502{
     1503  return 0;
    15161504}
    15171505void freeArrays1(ampl_info *)
     
    15211509{
    15221510}
    1523 void freeArgs(ampl_info * )
     1511void freeArgs(ampl_info *)
    15241512{
    15251513}
    15261514int ampl_obj_prec()
    15271515{
    1528     return 0;
    1529 }
    1530 void writeAmpl(ampl_info * )
     1516  return 0;
     1517}
     1518void writeAmpl(ampl_info *)
    15311519{
    15321520}
    15331521#endif
    15341522
     1523/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     1524*/
Note: See TracChangeset for help on using the changeset viewer.