source: trunk/Cbc/src/Cbc_ampl.cpp @ 775

Last change on this file since 775 was 775, checked in by ladanyi, 12 years ago

got rid of compiler warnings

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 43.5 KB
Line 
1/****************************************************************
2Copyright (C) 1997-2000 Lucent Technologies
3Modifications for Coin -  Copyright (C) 2006, International Business Machines Corporation and others.
4All Rights Reserved
5
6Permission to use, copy, modify, and distribute this software and
7its documentation for any purpose and without fee is hereby
8granted, provided that the above copyright notice appear in all
9copies and that both that the copyright notice and this
10permission notice and warranty disclaimer appear in supporting
11documentation, and that the name of Lucent or any of its entities
12not be used in advertising or publicity pertaining to
13distribution of the software without specific, written prior
14permission.
15
16LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
17INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
18IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
19SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
20WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
21IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
22ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
23THIS SOFTWARE.
24****************************************************************/
25#ifdef COIN_HAS_ASL
26
27#include "CbcConfig.h"
28#ifdef HAVE_UNISTD_H
29# include "unistd.h"
30#endif
31#include "CoinUtilsConfig.h"
32#include "CoinHelperFunctions.hpp"
33#include "CoinModel.hpp"
34#include "CoinSort.hpp"
35#include "CoinPackedMatrix.hpp"
36#include "CoinMpsIO.hpp"
37#include "CoinFloatEqual.hpp"
38#ifdef COIN_HAS_CLP
39#include "OsiClpSolverInterface.hpp"
40#endif
41#include "Cbc_ampl.h"
42extern "C" {
43# include "getstub.h"
44# include "asl_pfgh.h"
45}
46
47#include <string>
48#include <cassert>
49/* so decodePhrase and clpCheck can access */
50static ampl_info * saveInfo=NULL;
51// Set to 1 if algorithm found
52static char algFound[20]="";
53static char*
54checkPhrase(Option_Info *oi, keyword *kw, char *v)
55{
56  if (strlen(v))
57    printf("string %s\n",v);
58  // Say algorithm found
59  strcpy(algFound,kw->desc);
60  return v;
61}
62static char*
63checkPhrase2(Option_Info *oi, keyword *kw, char *v)
64{
65  if (strlen(v))
66    printf("string %s\n",v);
67  // put out keyword
68  saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+1)*sizeof(char *));
69  saveInfo->arguments[saveInfo->numberArguments++]=strdup(kw->desc);
70  return v;
71}
72static fint
73decodePhrase(char * phrase,ftnlen length)
74{
75  char * blank = strchr(phrase,' ');
76  if (blank) {
77    /* split arguments */
78    *blank='\0';
79    saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+2)*sizeof(char *));
80    saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
81    *blank=' ';
82    phrase=blank+1; /* move on */
83    if (strlen(phrase))
84      saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
85  } else if (strlen(phrase)) {
86    saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+1)*sizeof(char *));
87    saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
88  }
89  return 0;
90}
91static void
92sos_kludge(int nsos, int *sosbeg, double *sosref,int * sosind)
93{
94  // Adjust sosref if necessary to make monotonic increasing
95  int i, j, k;
96  // first sort
97  for (i=0;i<nsos;i++) {
98    k = sosbeg[i];
99    int end=sosbeg[i+1];
100    CoinSort_2(sosref+k,sosref+end,sosind+k);
101  }
102  double t, t1;
103  for(i = j = 0; i++ < nsos; ) {
104    k = sosbeg[i];
105    t = sosref[j];
106    while(++j < k) {
107      t1 = sosref[j];
108      t += 1e-10;
109      if (t1 <= t)
110        sosref[j] = t1 = t + 1e-10;
111      t = t1;
112    }
113  }
114}
115static char xxxxxx[20];
116#define VP (char*)
117static keyword keywds[] = { /* must be sorted */
118  { const_cast<char*>("barrier"),  checkPhrase,  (char *) xxxxxx ,
119    const_cast<char*>("-barrier")},
120  { const_cast<char*>("dual"),     checkPhrase,  (char *) xxxxxx ,
121    const_cast<char*>("-dualsimplex")},
122  { const_cast<char*>("help"),     checkPhrase2, (char *) xxxxxx ,
123    const_cast<char*>("-?")},
124  { const_cast<char*>("initial"),  checkPhrase,  (char *) xxxxxx ,
125    const_cast<char*>("-initialsolve")},
126  { const_cast<char*>("max"),      checkPhrase2, (char *) xxxxxx ,
127    const_cast<char*>("-maximize")},
128  { const_cast<char*>("maximize"), checkPhrase2, (char *) xxxxxx ,
129    const_cast<char*>("-maximize")},
130  { const_cast<char*>("primal"),   checkPhrase,  (char *) xxxxxx ,
131    const_cast<char*>("-primalsimplex")},
132  { const_cast<char*>("quit"),     checkPhrase2, (char *) xxxxxx ,
133    const_cast<char*>("-quit")},
134  { const_cast<char*>("wantsol"),  WS_val,       NULL,
135    const_cast<char*>("write .sol file (without -AMPL)")}
136};
137static Option_Info Oinfo = {
138  const_cast<char*>("cbc"),
139  const_cast<char*>("Cbc 1.04"),
140  const_cast<char*>("cbc_options"),
141  keywds,
142  nkeywds,
143  0,
144  const_cast<char*>(""),
145  0,
146  decodePhrase,
147  0,
148  0,
149  0,
150  20070606 };
151// strdup used to avoid g++ compiler warning
152static SufDecl suftab[] = {
153#if 0
154  { const_cast<char*>("current"), 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
155  { const_cast<char*>("current"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
156  { const_cast<char*>("direction"), 0, ASL_Sufkind_var },
157  { const_cast<char*>("down"), 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
158  { const_cast<char*>("down"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
159  { const_cast<char*>("priority"), 0, ASL_Sufkind_var },
160#endif
161  { const_cast<char*>("cut"), 0, ASL_Sufkind_con },
162  { const_cast<char*>("direction"), 0, ASL_Sufkind_var },
163  { const_cast<char*>("downPseudocost"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
164  { const_cast<char*>("priority"), 0, ASL_Sufkind_var },
165  { const_cast<char*>("ref"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
166  { const_cast<char*>("sos"), 0, ASL_Sufkind_var },
167  { const_cast<char*>("sos"), 0, ASL_Sufkind_con },
168  { const_cast<char*>("sosno"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
169  { const_cast<char*>("sosref"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
170  { const_cast<char*>("special"), 0, ASL_Sufkind_var },
171  { const_cast<char*>("special"), 0, ASL_Sufkind_con },
172  /*{ const_cast<char*>("special"), 0, ASL_Sufkind_con },*/
173  { const_cast<char*>("sstatus"), 0, ASL_Sufkind_var, 0 },
174  { const_cast<char*>("sstatus"), 0, ASL_Sufkind_con, 0 },
175  { const_cast<char*>("upPseudocost"), 0, ASL_Sufkind_var | ASL_Sufkind_real }
176#if 0
177  { const_cast<char*>("unbdd"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly},
178  { const_cast<char*>("up"), 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
179  { const_cast<char*>("up"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly }
180#endif
181};
182#include "float.h"
183#include "limits.h"
184static ASL *asl=NULL;
185static FILE *nl=NULL;
186
187static void
188mip_stuff(void)
189{
190  int i;
191  double *pseudoUp, *pseudoDown;
192  int *priority, *direction;
193  // To label cuts (there will be other uses for special)
194  int *cut;
195  // To label special variables - at present 1= must be >= 1 or <= -1
196  int * special;
197  SufDesc *dpup, *dpdown, *dpri, *ddir, *dcut, *dspecial;
198 
199  ddir = suf_get("direction", ASL_Sufkind_var);
200  direction = ddir->u.i;
201  dpri = suf_get("priority", ASL_Sufkind_var);
202  priority = dpri->u.i;
203  dspecial = suf_get("special", ASL_Sufkind_con);
204  dcut = suf_get("cut", ASL_Sufkind_con);
205  cut = dcut->u.i;
206  if (!cut) {
207    // try special
208    dcut = suf_get("special", ASL_Sufkind_con);
209    cut = dcut->u.i;
210  }
211  dspecial = suf_get("special", ASL_Sufkind_var);
212  special = dspecial->u.i;
213  dpdown = suf_get("downPseudocost", ASL_Sufkind_var);
214  pseudoDown = dpdown->u.r;
215  dpup = suf_get("upPseudocost", ASL_Sufkind_var);
216  pseudoUp = dpup->u.r;
217  assert(saveInfo);
218  int numberColumns = saveInfo->numberColumns;
219  if (direction) {
220    int baddir=0;
221    saveInfo->branchDirection = (int *) malloc(numberColumns*sizeof(int));
222    for (i=0;i<numberColumns;i++) {
223      int value = direction[i];
224      if (value<-1||value>1) {
225        baddir++;
226        value=0;
227      }
228      saveInfo->branchDirection[i]=value;
229    }
230    if (baddir)
231      fprintf(Stderr,
232              "Treating %d .direction values outside [-1, 1] as 0.\n",
233              baddir);
234  }
235  if (priority) {
236    int badpri=0;
237    saveInfo->priorities= (int *) malloc(numberColumns*sizeof(int));
238    for (i=0;i<numberColumns;i++) {
239      int value = priority[i];
240      if (value<0) {
241        badpri++;
242        value=0;
243      }
244      saveInfo->priorities[i]=value;
245    }
246    if (badpri)
247      fprintf(Stderr,
248              "Treating %d negative .priority values as 0\n",
249              badpri);
250  }
251  if (special) {
252    int badspecial=0;
253    saveInfo->special= (int *) malloc(numberColumns*sizeof(int));
254    for (i=0;i<numberColumns;i++) {
255      int value = special[i];
256      if (value<0) {
257        badspecial++;
258        value=0;
259      }
260      saveInfo->special[i]=value;
261    }
262    if (badspecial)
263      fprintf(Stderr,
264              "Treating %d negative special values as 0\n",
265              badspecial);
266  }
267  int numberRows = saveInfo->numberRows;
268  if (cut) {
269    int badcut=0;
270    saveInfo->cut= (int *) malloc(numberRows*sizeof(int));
271    for (i=0;i<numberRows;i++) {
272      int value = cut[i];
273      if (value<0) {
274        badcut++;
275        value=0;
276      }
277      saveInfo->cut[i]=value;
278    }
279    if (badcut)
280      fprintf(Stderr,
281              "Treating %d negative cut values as 0\n",
282              badcut);
283  }
284  if (pseudoDown||pseudoUp) {
285    int badpseudo=0;
286    if (!pseudoDown||!pseudoUp)
287      fprintf(Stderr,
288              "Only one set of pseudocosts - assumed same\n");
289    saveInfo->pseudoDown= (double *) malloc(numberColumns*sizeof(double));
290    saveInfo->pseudoUp = (double *) malloc(numberColumns*sizeof(double));
291    for (i=0;i<numberColumns;i++) {
292      double valueD=0.0, valueU=0.0;
293      if (pseudoDown) {
294        valueD = pseudoDown[i];
295        if (valueD<0) {
296          badpseudo++;
297          valueD=0.0;
298        }
299      }
300      if (pseudoUp) {
301        valueU = pseudoUp[i];
302        if (valueU<0) {
303          badpseudo++;
304          valueU=0.0;
305        }
306      }
307      if (!valueD)
308        valueD=valueU;
309      if (!valueU)
310        valueU=valueD;
311      saveInfo->pseudoDown[i]=valueD;
312      saveInfo->pseudoUp[i]=valueU;
313    }
314    if (badpseudo)
315      fprintf(Stderr,
316              "Treating %d negative pseudoCosts as 0.0\n",badpseudo);
317  }
318}
319static void
320stat_map(int *stat, int n, int *map, int mx, const char *what)
321{
322  int bad, i, i1=0, j, j1=0;
323  static char badfmt[] = "Coin driver: %s[%d] = %d\n";
324 
325  for(i = bad = 0; i < n; i++) {
326    if ((j = stat[i]) >= 0 && j <= mx)
327      stat[i] = map[j];
328    else {
329      stat[i] = 0;
330      i1 = i;
331      j1 = j;
332      if (!bad++)
333        fprintf(Stderr, badfmt, what, i, j);
334    }
335  }
336  if (bad > 1) {
337    if (bad == 2)
338      fprintf(Stderr, badfmt, what, i1, j1);
339    else
340      fprintf(Stderr,
341              "Coin driver: %d messages about bad %s values suppressed.\n",
342              bad-1, what);
343  }
344}
345int
346readAmpl(ampl_info * info, int argc, char **argv, void ** coinModel)
347{
348  char *stub;
349  ograd *og;
350  int i;
351  SufDesc *csd;
352  SufDesc *rsd;
353  /*bool *basis, *lower;*/
354  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
355  char * environment = getenv("cbc_options");
356  char tempBuffer[20];
357  double * obj;
358  double * columnLower;
359  double * columnUpper;
360  double * rowLower;
361  double * rowUpper;
362  char ** saveArgv=argv;
363  char fileName[1000];
364  if (argc>1)
365    strcpy(fileName,argv[1]);
366  else
367    fileName[0]='\0';
368  int nonLinearType=-1;
369  // testosi parameter - if >= 10 then go in through coinModel
370  for (i=1;i<argc;i++) {
371    if (!strncmp(argv[i],"testosi",7)) {
372      char * equals = strchr(argv[i],'=');
373      if (equals&&atoi(equals+1)>=10&&atoi(equals+1)<=20) {
374        nonLinearType=atoi(equals+1);
375        break;
376      }
377    }
378  }
379  int saveArgc = argc;
380  if (info->numberRows!=-1234567)
381    memset(info,0,sizeof(ampl_info)); // overwrite unless magic number set
382  /* save so can be accessed by decodePhrase */
383  saveInfo = info;
384  info->numberArguments=0;
385  info->arguments=(char **) malloc(2*sizeof(char *));
386  info->arguments[info->numberArguments++]=strdup("ampl");
387  info->arguments[info->numberArguments++]=strdup("cbc");
388  asl = ASL_alloc(ASL_read_f);
389  stub = getstub(&argv, &Oinfo);
390  if (!stub)
391    usage_ASL(&Oinfo, 1);
392  nl = jac0dim(stub, 0);
393  suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
394 
395  /* set A_vals to get the constraints column-wise (malloc so can be freed) */
396  A_vals = (double *) malloc(nzc*sizeof(double));
397  if (!A_vals) {
398    printf("no memory\n");
399    return 1;
400  }
401  /* say we want primal solution */
402  want_xpi0=1;
403  /* for basis info */
404  info->columnStatus = (int *) malloc(n_var*sizeof(int));
405  info->rowStatus = (int *) malloc(n_con*sizeof(int));
406  csd = suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
407  rsd = suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
408  if (!(nlvc+nlvo)&&nonLinearType<10) {
409    /* read linear model*/
410    f_read(nl,0);
411    // see if any sos
412    if (true) {
413      char *sostype;
414      int nsosnz, *sosbeg, *sosind, * sospri;
415      double *sosref;
416      int nsos;
417      int i = ASL_suf_sos_explict_free;
418      int copri[2], **p_sospri;
419      copri[0] = 0;
420      copri[1] = 0;
421      p_sospri = &sospri;
422      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
423                     &sosbeg, &sosind, &sosref);
424      if (nsos) {
425        info->numberSos=nsos;
426        info->sosType = (char *) malloc(nsos);
427        info->sosPriority = (int *) malloc(nsos*sizeof(int));
428        info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
429        info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
430        info->sosReference = (double *) malloc(nsosnz*sizeof(double));
431        sos_kludge(nsos, sosbeg, sosref,sosind);
432        for (int i=0;i<nsos;i++) {
433          int ichar = sostype[i];
434          assert (ichar=='1'||ichar=='2');
435          info->sosType[i]=ichar-'0';
436        }       
437        memcpy(info->sosPriority,sospri,nsos*sizeof(int));
438        memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
439        memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
440        memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
441      }
442    }
443   
444    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
445    Oinfo.uinfo = tempBuffer;
446    if (getopts(argv, &Oinfo))
447      return 1;
448    /* objective*/
449    obj = (double *) malloc(n_var*sizeof(double));
450    for (i=0;i<n_var;i++)
451      obj[i]=0.0;
452    if (n_obj) {
453      for (og = Ograd[0];og;og = og->next)
454        obj[og->varno] = og->coef;
455    }
456    if (objtype[0])
457      info->direction=-1.0;
458    else
459      info->direction=1.0;
460    info->offset=objconst(0);
461    /* Column bounds*/
462    columnLower = (double *) malloc(n_var*sizeof(double));
463    columnUpper = (double *) malloc(n_var*sizeof(double));
464#define COIN_DBL_MAX DBL_MAX
465    for (i=0;i<n_var;i++) {
466      columnLower[i]=LUv[2*i];
467      if (columnLower[i]<= negInfinity)
468        columnLower[i]=-COIN_DBL_MAX;
469      columnUpper[i]=LUv[2*i+1];
470      if (columnUpper[i]>= Infinity)
471        columnUpper[i]=COIN_DBL_MAX;
472    }
473    /* Row bounds*/
474    rowLower = (double *) malloc(n_con*sizeof(double));
475    rowUpper = (double *) malloc(n_con*sizeof(double));
476    for (i=0;i<n_con;i++) {
477      rowLower[i]=LUrhs[2*i];
478      if (rowLower[i]<= negInfinity)
479        rowLower[i]=-COIN_DBL_MAX;
480      rowUpper[i]=LUrhs[2*i+1];
481      if (rowUpper[i]>= Infinity)
482        rowUpper[i]=COIN_DBL_MAX;
483    }
484    info->numberRows=n_con;
485    info->numberColumns=n_var;
486    info->numberElements=nzc;
487    info->numberBinary=nbv;
488    info->numberIntegers=niv+nbv;
489    info->objective=obj;
490    info->rowLower=rowLower;
491    info->rowUpper=rowUpper;
492    info->columnLower=columnLower;
493    info->columnUpper=columnUpper;
494    info->starts=A_colstarts;
495    /*A_colstarts=NULL;*/
496    info->rows=A_rownos;
497    /*A_rownos=NULL;*/
498    info->elements=A_vals;
499    /*A_vals=NULL;*/
500    info->primalSolution=NULL;
501    /* put in primalSolution if exists */
502    if (X0) {
503      info->primalSolution=(double *) malloc(n_var*sizeof(double));
504      memcpy(info->primalSolution,X0,n_var*sizeof(double));
505    }
506    info->dualSolution=NULL;
507    if (niv+nbv>0)
508      mip_stuff(); // get any extra info
509    if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
510        ||(rsd->kind & ASL_Sufkind_input)) {
511      /* convert status - need info on map */
512      static int map[] = {1, 3, 1, 1, 2, 1, 1};
513      stat_map(info->columnStatus, n_var, map, 6, "incoming columnStatus");
514      stat_map(info->rowStatus, n_con, map, 6, "incoming rowStatus");
515    } else {
516      /* all slack basis */
517      // leave status for output */
518#if 0
519      free(info->rowStatus);
520      info->rowStatus=NULL;
521      free(info->columnStatus);
522      info->columnStatus=NULL;
523#endif
524    }
525  } else {
526    // QP
527    // Add .nl if not there
528    if (!strstr(fileName,".nl"))
529      strcat(fileName,".nl");
530    CoinModel * model = new CoinModel((nonLinearType>10) ? 2 : 1,fileName,info);
531    if (model->numberRows()>0||model->numberColumns()>0)
532      *coinModel=(void *) model;
533    Oinfo.uinfo = tempBuffer;
534    if (getopts(argv, &Oinfo))
535      return 1;
536    Oinfo.wantsol=1;
537    if (objtype[0])
538      info->direction=-1.0;
539    else
540      info->direction=1.0;
541    model->setOptimizationDirection(info->direction);
542    info->offset=objconst(0);
543    info->numberRows=n_con;
544    info->numberColumns=n_var;
545    info->numberElements=nzc;
546    info->numberBinary=nbv;
547    int numberIntegers=niv+nlvci+nlvoi+nbv;
548    if (nlvci+nlvoi+nlvc+nlvo) {
549      // Non linear
550      // No idea if there are overlaps so compute
551      int numberIntegers=0;
552      for ( i=0; i<n_var;i++) {
553        if (model->columnIsInteger(i))
554          numberIntegers++;
555      }
556    }
557    info->numberIntegers=numberIntegers;
558    // Say nonlinear if it is
559    info->nonLinear=nlvc+nlvo;
560    if (numberIntegers>0) {
561      mip_stuff(); // get any extra info
562      if (info->cut) 
563        model->setCutMarker(info->numberRows,info->cut);
564      if (info->priorities)
565        model->setPriorities(info->numberColumns,info->priorities);
566    }
567  }
568  /* add -solve - unless something there already
569   - also check for sleep=yes */
570  {
571    int found=0;
572    int foundLog=0;
573    int foundSleep=0;
574    const char * something[]={"solve","branch","duals","primals","user"};
575    for (i=0;i<info->numberArguments;i++) {
576      unsigned int j;
577      const char * argument = info->arguments[i];
578      for (j=0;j<sizeof(something)/sizeof(char *);j++) {
579        const char * check = something[j];
580        if (!strncmp(argument,check,sizeof(check))) {
581          found=(int)(j+1);
582        } else if (!strncmp(argument,"log",3)) {
583          foundLog=1;
584        } else if (!strncmp(argument,"sleep",5)) {
585          foundSleep=1;
586        }
587      }
588    }
589    if (foundLog) {
590      /* print options etc */
591      for (i=0;i<saveArgc;i++)
592        printf("%s ",saveArgv[i]);
593      printf("\n");
594      if (environment)
595        printf("env %s\n",environment);
596      /*printf("%d rows %d columns %d elements\n",n_con,n_var,nzc);*/
597    }
598    if (!found) {
599      if (!strlen(algFound)) {
600        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
601        info->arguments[info->numberArguments++]=strdup("-solve");
602      } else {
603        // use algorithm from keyword
604        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
605        info->arguments[info->numberArguments++]=strdup(algFound);
606      }
607    }
608    if (foundSleep) {
609      /* let user copy .nl file */
610      fprintf(stderr,"You can copy .nl file %s for debug purposes or attach debugger\n",saveArgv[1]);
611      fprintf(stderr,"Type q to quit, anything else to continue\n");
612      char getChar = getc(stdin);
613      if (getChar=='q'||getChar=='Q')
614        exit(1);
615    }
616  }
617  /* add -quit */
618  info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
619  info->arguments[info->numberArguments++]=strdup("-quit");
620  return 0;
621}
622void freeArrays1(ampl_info * info)
623{
624  free(info->objective);
625  info->objective=NULL;
626  free(info->rowLower);
627  info->rowLower=NULL;
628  free(info->rowUpper);
629  info->rowUpper=NULL;
630  free(info->columnLower);
631  info->columnLower=NULL;
632  free(info->columnUpper);
633  info->columnUpper=NULL;
634  /* this one not freed by ASL_free */
635  free(info->elements);
636  info->elements=NULL;
637  free(info->primalSolution);
638  info->primalSolution=NULL;
639  free(info->dualSolution);
640  info->dualSolution=NULL;
641  /*free(info->rowStatus);
642  info->rowStatus=NULL;
643  free(info->columnStatus);
644  info->columnStatus=NULL;*/
645}
646void freeArrays2(ampl_info * info)
647{
648  free(info->primalSolution);
649  info->primalSolution=NULL;
650  free(info->dualSolution);
651  info->dualSolution=NULL;
652  free(info->rowStatus);
653  info->rowStatus=NULL;
654  free(info->columnStatus);
655  info->columnStatus=NULL;
656  free(info->priorities);
657  info->priorities=NULL;
658  free(info->branchDirection);
659  info->branchDirection=NULL;
660  free(info->pseudoDown);
661  info->pseudoDown=NULL;
662  free(info->pseudoUp);
663  info->pseudoUp=NULL;
664  free(info->sosType);
665  info->sosType=NULL;
666  free(info->sosPriority);
667  info->sosPriority=NULL;
668  free(info->sosStart);
669  info->sosStart=NULL;
670  free(info->sosIndices);
671  info->sosIndices=NULL;
672  free(info->sosReference);
673  info->sosReference=NULL;
674  free(info->cut);
675  info->cut=NULL;
676  ASL_free(&asl);
677}
678void freeArgs(ampl_info * info)
679{
680  int i;
681  for ( i=0; i<info->numberArguments;i++)
682    free(info->arguments[i]);
683  free(info->arguments);
684}
685int ampl_obj_prec()
686{
687  return obj_prec();
688}
689void writeAmpl(ampl_info * info)
690{
691  char buf[1000];
692  typedef struct { const char *msg; int code; int wantObj; } Sol_info;
693  static Sol_info solinfo[] = {
694    { "optimal solution",                       000, 1 },
695    { "infeasible",                             200, 1 },
696    { "unbounded",                              300, 0 },
697    { "iteration limit etc",                    400, 1 },
698    { "solution limit",                         401, 1 },
699    { "ran out of space",                       500, 0 },
700    { "status unknown",                         501, 1 },
701    { "bug!",                                   502, 0 },
702    { "best MIP solution so far restored",      101, 1 },
703    { "failed to restore best MIP solution",    503, 1 },
704    { "optimal (?) solution",                   100, 1 }
705  };
706  /* convert status - need info on map */
707  static int map[] = {0, 3, 4, 1};
708  sprintf(buf,"%s %s",Oinfo.bsname,info->buffer);
709  solve_result_num = solinfo[info->problemStatus].code;
710  if (info->columnStatus) {
711    stat_map(info->columnStatus, n_var, map, 4, "outgoing columnStatus");
712    stat_map(info->rowStatus, n_con, map, 4, "outgoing rowStatus");
713    suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
714    suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
715  }
716  write_sol(buf,info->primalSolution,info->dualSolution,&Oinfo);
717}
718/* Read a problem from AMPL nl file
719 */
720CoinModel::CoinModel( int nonLinear, const char * fileName,const void * info)
721 :  numberRows_(0),
722    maximumRows_(0),
723    numberColumns_(0),
724    maximumColumns_(0),
725    numberElements_(0),
726    maximumElements_(0),
727    numberQuadraticElements_(0),
728    maximumQuadraticElements_(0),
729    optimizationDirection_(1.0),
730    objectiveOffset_(0.0),
731    rowLower_(NULL),
732    rowUpper_(NULL),
733    rowType_(NULL),
734    objective_(NULL),
735    columnLower_(NULL),
736    columnUpper_(NULL),
737    integerType_(NULL),
738    columnType_(NULL),
739    start_(NULL),
740    elements_(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    logLevel_(0),
756    type_(-1),
757    links_(0)
758{
759  problemName_ = strdup("");
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  free(problemName_);
1270  problemName_=strdup("???");
1271
1272  // Build by row from scratch
1273  const double * element = matrixByRow.getElements();
1274  const int * column = matrixByRow.getIndices();
1275  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
1276  const int * rowLength = matrixByRow.getVectorLengths();
1277  for (i=0;i<numberRows;i++) {
1278    addRow(rowLength[i],column+rowStart[i],
1279           element+rowStart[i],rowLower[i],rowUpper[i]);
1280  }
1281  // Now do column part
1282  for (i=0;i<numberColumns;i++) {
1283    setColumnBounds(i,columnLower[i],columnUpper[i]);
1284    setColumnObjective(i,objective[i]);
1285  }
1286  for ( i=numberColumns-numberBinary-numberIntegers;
1287        i<numberColumns;i++) {
1288    setColumnIsInteger(i,true);
1289  }
1290  // and non linear
1291  for (i=numberAllNonLinearBoth-numberIntegerNonLinearBoth;
1292       i<numberAllNonLinearBoth;i++) {
1293    setColumnIsInteger(i,true);
1294  }
1295  for (i=numberAllNonLinearConstraints-numberIntegerNonLinearConstraints;
1296       i<numberAllNonLinearConstraints;i++) {
1297    setColumnIsInteger(i,true);
1298  }
1299  for (i=numberAllNonLinearObjective-numberIntegerNonLinearObjective;
1300       i<numberAllNonLinearObjective;i++) {
1301    setColumnIsInteger(i,true);
1302  }
1303  free(columnLower);
1304  free(columnUpper);
1305  free(rowLower);
1306  free(rowUpper);
1307  free(objective);
1308  // do names
1309  int iRow;
1310  for (iRow=0;iRow<numberRows_;iRow++) {
1311    char name[9];
1312    sprintf(name,"r%7.7d",iRow);
1313    setRowName(iRow,name);
1314  }
1315  int iColumn;
1316  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1317    char name[9];
1318    sprintf(name,"c%7.7d",iColumn);
1319    setColumnName(iColumn,name);
1320  }
1321  if (colqp) {
1322    // add in quadratic
1323    int nz = 1 + n_con;
1324    int nOdd=0;
1325    fint ** rowqp = colqp + nz;
1326    double ** delsqp = (double **)(rowqp + nz);
1327    for (i=0;i<=n_con;i++) {
1328      int nels = z[i];
1329      if (nels) {
1330        double * element = delsqp[i];
1331        int * start = (int *) colqp[i];
1332        int * row = (int *) rowqp[i];
1333        if (!element) {
1334          // odd row - probably not quadratic
1335          nOdd++;
1336          continue;
1337        }
1338#if 0
1339        printf("%d quadratic els\n",nels);
1340        for (int j=0;j<n_var;j++) {
1341          for (int k=start[j];k<start[j+1];k++)
1342            printf("%d %d %g\n",j,row[k],element[k]);
1343        }
1344#endif
1345        if (i) {
1346          int iRow = i-1;
1347          for (int j=0;j<n_var;j++) {
1348            for (int k=start[j];k<start[j+1];k++) {
1349              int kColumn = row[k];
1350              double value = element[k];
1351              // ampl gives twice with assumed 0.5
1352              if (kColumn<j)
1353                continue;
1354              else if (kColumn==j)
1355                value *= 0.5;
1356              const char * expr = getElementAsString(iRow,j);
1357              double constant=0.0;
1358              bool linear;
1359              if (expr&&strcmp(expr,"Numeric")) {
1360                linear=false;
1361              } else {
1362                constant = getElement(iRow,j);
1363                linear=true;
1364              }
1365              char temp[1000];
1366              char temp2[30];
1367              if (value==1.0) 
1368                sprintf(temp2,"c%7.7d",kColumn);
1369              else
1370                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1371              if (linear) {
1372                if (!constant)
1373                  strcpy(temp,temp2);
1374                else if (value>0.0) 
1375                  sprintf(temp,"%g+%s",constant,temp2);
1376                else
1377                  sprintf(temp,"%g%s",constant,temp2);
1378              } else {
1379                if (value>0.0) 
1380                  sprintf(temp,"%s+%s",expr,temp2);
1381                else
1382                  sprintf(temp,"%s%s",expr,temp2);
1383              }
1384              assert (strlen(temp)<1000);
1385              setElement(iRow,j,temp);
1386              if (amplInfo->logLevel>1)
1387                printf("el for row %d column c%7.7d is %s\n",iRow,j,temp);
1388            }
1389          }
1390        } else {
1391          // objective
1392          for (int j=0;j<n_var;j++) {
1393            for (int k=start[j];k<start[j+1];k++) {
1394              int kColumn = row[k];
1395              double value = element[k];
1396              // ampl gives twice with assumed 0.5
1397              if (kColumn<j)
1398                continue;
1399              else if (kColumn==j)
1400                value *= 0.5;
1401              const char * expr = getColumnObjectiveAsString(j);
1402              double constant=0.0;
1403              bool linear;
1404              if (expr&&strcmp(expr,"Numeric")) {
1405                linear=false;
1406              } else {
1407                constant = getColumnObjective(j);
1408                linear=true;
1409              }
1410              char temp[1000];
1411              char temp2[30];
1412              if (value==1.0) 
1413                sprintf(temp2,"c%7.7d",kColumn);
1414              else
1415                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1416              if (linear) {
1417                if (!constant)
1418                  strcpy(temp,temp2);
1419                else if (value>0.0) 
1420                  sprintf(temp,"%g+%s",constant,temp2);
1421                else
1422                  sprintf(temp,"%g%s",constant,temp2);
1423              } else {
1424                if (value>0.0) 
1425                  sprintf(temp,"%s+%s",expr,temp2);
1426                else
1427                  sprintf(temp,"%s%s",expr,temp2);
1428              }
1429              assert (strlen(temp)<1000);
1430              setObjective(j,temp);
1431              if (amplInfo->logLevel>1)
1432                printf("el for objective column c%7.7d is %s\n",j,temp);
1433            }
1434          }
1435        }
1436      }
1437    }
1438    if (nOdd) {
1439      printf("%d non-linear constraints could not be converted to quadratic\n",nOdd);
1440      exit(77);
1441    }
1442  }
1443  free(colqp);
1444  free(z);
1445  // see if any sos
1446  {
1447    char *sostype;
1448    int nsosnz, *sosbeg, *sosind, * sospri;
1449    double *sosref;
1450    int nsos;
1451    int i = ASL_suf_sos_explict_free;
1452    int copri[2], **p_sospri;
1453    copri[0] = 0;
1454    copri[1] = 0;
1455    p_sospri = &sospri;
1456    nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1457                   &sosbeg, &sosind, &sosref);
1458    if (nsos) {
1459      numberSOS_=nsos;
1460      typeSOS_ = new int [numberSOS_];
1461      prioritySOS_ = new int [numberSOS_];
1462      startSOS_ = new int [numberSOS_+1];
1463      memberSOS_ = new int[nsosnz];
1464      referenceSOS_ = new double [nsosnz];
1465      sos_kludge(nsos, sosbeg, sosref,sosind);
1466      for (int i=0;i<nsos;i++) {
1467        int ichar = sostype[i];
1468        assert (ichar=='1'||ichar=='2');
1469        typeSOS_[i]=ichar-'0';
1470      } 
1471      memcpy(prioritySOS_,sospri,nsos*sizeof(int));
1472      memcpy(startSOS_,sosbeg,(nsos+1)*sizeof(int));
1473      memcpy(memberSOS_,sosind,nsosnz*sizeof(int));
1474      memcpy(referenceSOS_,sosref,nsosnz*sizeof(double));
1475    }
1476  }
1477}
1478#else
1479#include "Cbc_ampl.h"
1480int
1481readAmpl(ampl_info * info, int argc, char **argv, void ** coinModel)
1482{
1483  return 0;
1484}
1485void freeArrays1(ampl_info * info)
1486{
1487}
1488void freeArrays2(ampl_info * info)
1489{
1490}
1491void freeArgs(ampl_info * info)
1492{
1493}
1494int ampl_obj_prec()
1495{
1496  return 0;
1497}
1498void writeAmpl(ampl_info * info)
1499{
1500}
1501#endif
Note: See TracBrowser for help on using the repository browser.