source: trunk/Clp/src/Clp_ampl.cpp @ 2470

Last change on this file since 2470 was 2385, checked in by unxusr, 8 months ago

formatting

File size: 46.7 KB
Line 
1/* $Id: Clp_ampl.cpp 2030 2014-04-14 17:34:30Z forrest $ */
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 Clp_ampl.cpp
28
29  Interface routines for AMPL.
30*/
31
32#include "ClpConfig.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#include "Clp_ampl.h"
47extern "C" {
48#include "getstub.h"
49#include "asl_pfgh.h"
50}
51
52#include <string>
53#include <cassert>
54/* so decodePhrase and clpCheck can access */
55static ampl_info *saveInfo = NULL;
56// Set to 1 if algorithm found
57static char algFound[20] = "";
58static char *
59checkPhrase(Option_Info *oi, keyword *kw, char *v)
60{
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 *
68checkPhrase2(Option_Info *oi, keyword *kw, char *v)
69{
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;
76}
77static fint
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;
95}
96static void
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  }
119}
120static char xxxxxx[20];
121#define VP (char *)
122static 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)") }
141};
142static 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
156};
157// strdup used to avoid g++ compiler warning
158static SufDecl suftab[] = {
159#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 },
166#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 }
182#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 }
186#endif
187};
188#include "float.h"
189#include "limits.h"
190static ASL *asl = NULL;
191static FILE *nl = NULL;
192
193static void
194mip_stuff(void)
195{
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);
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;
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;
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  }
324}
325static void
326stat_map(int *stat, int n, int *map, int mx, const char *what)
327{
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');
450        }
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;
472    else
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;
603        }
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
871    asl = ASL_alloc(ASL_read_f);
872    nl = jac0dim(stub, 0);
873    free(stub);
874    suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
875
876    /* set A_vals to get the constraints column-wise (malloc so can be freed) */
877    A_vals = (double *)malloc(nzc * sizeof(double));
878    if (!A_vals) {
879      printf("no memory\n");
880      return;
881    }
882    /* say we want primal solution */
883    want_xpi0 = 1;
884    /* for basis info */
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);
892#else
893    f_read(nl, 0);
894#endif
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);
1409            }
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          }
1456        }
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  }
1498}
1499#else
1500#include "Clp_ampl.h"
1501int readAmpl(ampl_info *, int, char **, void **)
1502{
1503  return 0;
1504}
1505void freeArrays1(ampl_info *)
1506{
1507}
1508void freeArrays2(ampl_info *)
1509{
1510}
1511void freeArgs(ampl_info *)
1512{
1513}
1514int ampl_obj_prec()
1515{
1516  return 0;
1517}
1518void writeAmpl(ampl_info *)
1519{
1520}
1521#endif
1522
1523/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
1524*/
Note: See TracBrowser for help on using the repository browser.