source: trunk/Cbc/src/Cbc_ampl.cpp

Last change on this file was 2616, checked in by stefan, 5 months ago

fix argument processing in ampl interface

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