source: branches/sandbox/Cbc/src/Cbc_ampl.cpp @ 1275

Last change on this file since 1275 was 1275, checked in by bjarni, 10 years ago

Test commit from Bjarni

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