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

Last change on this file since 789 was 789, checked in by forrest, 12 years ago

trying to move ampl stuff away

  • 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    moreInfo_(NULL),
756    logLevel_(0),
757    type_(-1),
758    links_(0)
759{
760  problemName_ = strdup("");
761  int status=0;
762  if (!strcmp(fileName,"-")||!strcmp(fileName,"stdin")) {
763    // stdin
764  } else {
765    std::string name=fileName;
766    bool readable = fileCoinReadable(name);
767    if (!readable) {
768      std::cerr<<"Unable to open file "
769               <<fileName<<std::endl;
770      status = -1;
771    }
772  }
773  if (!status) {
774    gdb(nonLinear,fileName,info);
775  }
776}
777#if 0
778static real
779qterm(ASL *asl, fint *colq, fint *rowq, real *delsq)
780{
781  double t, t1, *x, *x0, *xe;
782  fint *rq0, *rqe;
783 
784  t = 0.;
785  x0 = x = X0;
786  xe = x + n_var;
787  rq0 = rowq;
788  while(x < xe) {
789    t1 = *x++;
790    rqe = rq0 + *++colq;
791    while(rowq < rqe)
792      t += t1*x0[*rowq++]**delsq++;
793  }
794  return 0.5 * t;
795}
796#endif
797// stolen from IPopt with changes
798typedef struct {
799  double obj_sign_;
800  ASL_pfgh * asl_;
801  double * non_const_x_;
802  int * column_; // for jacobian
803  int * rowStart_;
804  double * gradient_;
805  double * constraintValues_;
806  int nz_h_full_; // number of nonzeros in hessian
807  int nerror_;
808  bool objval_called_with_current_x_;
809  bool conval_called_with_current_x_;
810  bool jacval_called_with_current_x_;
811} CbcAmplInfo;
812
813void
814CoinModel::gdb( int nonLinear, const char * fileName,const void * info)
815{
816  const ampl_info * amplInfo = (const ampl_info *) info;
817  ograd *og=NULL;
818  int i;
819  SufDesc *csd=NULL;
820  SufDesc *rsd=NULL;
821  /*bool *basis, *lower;*/
822  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
823  //char tempBuffer[20];
824  double * objective=NULL;
825  double * columnLower=NULL;
826  double * columnUpper=NULL;
827  double * rowLower=NULL;
828  double * rowUpper=NULL;
829  int * columnStatus=NULL;
830  int * rowStatus=NULL;
831  int numberRows=-1;
832  int numberColumns=-1;
833  int numberElements=-1;
834  int numberBinary=-1;
835  int numberIntegers=-1;
836  int numberAllNonLinearBoth=0;
837  int numberIntegerNonLinearBoth=0;
838  int numberAllNonLinearConstraints=0;
839  int numberIntegerNonLinearConstraints=0;
840  int numberAllNonLinearObjective=0;
841  int numberIntegerNonLinearObjective=0;
842  double * primalSolution=NULL;
843  double direction=1.0;
844  char * stub = strdup(fileName);
845  CoinPackedMatrix matrixByRow;
846  fint ** colqp=NULL;
847  int *z = NULL;
848  if (nonLinear==0) {
849    // linear
850    asl = ASL_alloc(ASL_read_f);
851    nl = jac0dim(stub, 0);
852    free(stub);
853    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
854   
855    /* set A_vals to get the constraints column-wise (malloc so can be freed) */
856    A_vals = (double *) malloc(nzc*sizeof(double));
857    if (!A_vals) {
858      printf("no memory\n");
859      return ;
860    }
861    /* say we want primal solution */
862    want_xpi0=1;
863    /* for basis info */
864    columnStatus = (int *) malloc(n_var*sizeof(int));
865    rowStatus = (int *) malloc(n_con*sizeof(int));
866    csd = suf_iput("sstatus", ASL_Sufkind_var, columnStatus);
867    rsd = suf_iput("sstatus", ASL_Sufkind_con, rowStatus);
868    /* read linear model*/
869    f_read(nl,0);
870    // see if any sos
871    if (true) {
872      char *sostype;
873      int nsosnz, *sosbeg, *sosind, * sospri;
874      double *sosref;
875      int nsos;
876      int i = ASL_suf_sos_explict_free;
877      int copri[2], **p_sospri;
878      copri[0] = 0;
879      copri[1] = 0;
880      p_sospri = &sospri;
881      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
882                     &sosbeg, &sosind, &sosref);
883      if (nsos) {
884        abort();
885#if 0
886        info->numberSos=nsos;
887        info->sosType = (char *) malloc(nsos);
888        info->sosPriority = (int *) malloc(nsos*sizeof(int));
889        info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
890        info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
891        info->sosReference = (double *) malloc(nsosnz*sizeof(double));
892        sos_kludge(nsos, sosbeg, sosref,sosind);
893        for (int i=0;i<nsos;i++) {
894          int ichar = sostype[i];
895          assert (ichar=='1'||ichar=='2');
896          info->sosType[i]=ichar-'0';
897        }       
898        memcpy(info->sosPriority,sospri,nsos*sizeof(int));
899        memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
900        memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
901        memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
902#endif
903      }
904    }
905   
906    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
907    //Oinfo.uinfo = tempBuffer;
908    //if (getopts(argv, &Oinfo))
909    //return 1;
910    /* objective*/
911    objective = (double *) malloc(n_var*sizeof(double));
912    for (i=0;i<n_var;i++)
913      objective[i]=0.0;
914    if (n_obj) {
915      for (og = Ograd[0];og;og = og->next)
916        objective[og->varno] = og->coef;
917    }
918    if (objtype[0])
919      direction=-1.0;
920    else
921      direction=1.0;
922    objectiveOffset_=objconst(0);
923    /* Column bounds*/
924    columnLower = (double *) malloc(n_var*sizeof(double));
925    columnUpper = (double *) malloc(n_var*sizeof(double));
926#define COIN_DBL_MAX DBL_MAX
927    for (i=0;i<n_var;i++) {
928      columnLower[i]=LUv[2*i];
929      if (columnLower[i]<= negInfinity)
930        columnLower[i]=-COIN_DBL_MAX;
931      columnUpper[i]=LUv[2*i+1];
932      if (columnUpper[i]>= Infinity)
933        columnUpper[i]=COIN_DBL_MAX;
934    }
935    /* Row bounds*/
936    rowLower = (double *) malloc(n_con*sizeof(double));
937    rowUpper = (double *) malloc(n_con*sizeof(double));
938    for (i=0;i<n_con;i++) {
939      rowLower[i]=LUrhs[2*i];
940      if (rowLower[i]<= negInfinity)
941        rowLower[i]=-COIN_DBL_MAX;
942      rowUpper[i]=LUrhs[2*i+1];
943      if (rowUpper[i]>= Infinity)
944        rowUpper[i]=COIN_DBL_MAX;
945    }
946    numberRows=n_con;
947    numberColumns=n_var;
948    numberElements=nzc;
949    numberBinary=nbv;
950    numberIntegers=niv;
951    /* put in primalSolution if exists */
952    if (X0) {
953      primalSolution=(double *) malloc(n_var*sizeof(double));
954      memcpy( primalSolution,X0,n_var*sizeof(double));
955    }
956    //double * dualSolution=NULL;
957    if (niv+nbv>0)
958      mip_stuff(); // get any extra info
959    if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
960        ||(rsd->kind & ASL_Sufkind_input)) {
961      /* convert status - need info on map */
962      static int map[] = {1, 3, 1, 1, 2, 1, 1};
963      stat_map(columnStatus, n_var, map, 6, "incoming columnStatus");
964      stat_map(rowStatus, n_con, map, 6, "incoming rowStatus");
965    } else {
966      /* all slack basis */
967      // leave status for output */
968#if 0
969      free(rowStatus);
970      rowStatus=NULL;
971      free(columnStatus);
972      columnStatus=NULL;
973#endif
974    }
975    CoinPackedMatrix columnCopy(true,numberRows,numberColumns,numberElements,
976                                A_vals,A_rownos,A_colstarts,NULL);
977    matrixByRow.reverseOrderedCopyOf(columnCopy);
978  } else if (nonLinear==1) {
979    // quadratic
980    asl = ASL_alloc(ASL_read_fg);
981    nl = jac0dim(stub, (long) strlen(stub));
982    free(stub);
983    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
984    /* read  model*/
985    X0 = (double*) malloc(n_var*sizeof(double));
986    CoinZeroN(X0,n_var);
987    qp_read(nl,0);
988    assert (n_obj==1);
989    int nz = 1 + n_con;
990    colqp = (fint**) malloc(nz*(2*sizeof(int*)
991                                      + sizeof(double*)));
992    fint ** rowqp = colqp + nz;
993    double ** delsqp = (double **)(rowqp + nz);
994    z = (int*) malloc(nz*sizeof(int));
995    for (i=0;i<=n_con;i++) {
996      z[i] = nqpcheck(-i, rowqp+i, colqp+i, delsqp+i);
997    }
998    qp_opify();
999    /* objective*/
1000    objective = (double *) malloc(n_var*sizeof(double));
1001    for (i=0;i<n_var;i++)
1002      objective[i]=0.0;
1003    if (n_obj) {
1004      for (og = Ograd[0];og;og = og->next)
1005        objective[og->varno] = og->coef;
1006    }
1007    if (objtype[0])
1008      direction=-1.0;
1009    else
1010      direction=1.0;
1011    objectiveOffset_=objconst(0);
1012    /* Column bounds*/
1013    columnLower = (double *) malloc(n_var*sizeof(double));
1014    columnUpper = (double *) malloc(n_var*sizeof(double));
1015    for (i=0;i<n_var;i++) {
1016      columnLower[i]=LUv[2*i];
1017      if (columnLower[i]<= negInfinity)
1018        columnLower[i]=-COIN_DBL_MAX;
1019      columnUpper[i]=LUv[2*i+1];
1020      if (columnUpper[i]>= Infinity)
1021        columnUpper[i]=COIN_DBL_MAX;
1022    }
1023    // Build by row from scratch
1024    //matrixByRow.reserve(n_var,nzc,true);
1025    // say row orderded
1026    matrixByRow.transpose();
1027    /* Row bounds*/
1028    rowLower = (double *) malloc(n_con*sizeof(double));
1029    rowUpper = (double *) malloc(n_con*sizeof(double));
1030    CoinBigIndex * rowStart = new CoinBigIndex [n_con+1];
1031    int * column = new int [nzc];
1032    double * element = new double [nzc];
1033    rowStart[0]=0;
1034    numberElements=0;
1035    for (i=0;i<n_con;i++) {
1036      rowLower[i]=LUrhs[2*i];
1037      if (rowLower[i]<= negInfinity)
1038        rowLower[i]=-COIN_DBL_MAX;
1039      rowUpper[i]=LUrhs[2*i+1];
1040      if (rowUpper[i]>= Infinity)
1041        rowUpper[i]=COIN_DBL_MAX;
1042      for(cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
1043        column[numberElements]=cg->varno;
1044        element[numberElements++]=cg->coef;
1045      }
1046      rowStart[i+1]=numberElements;
1047    }
1048    assert (numberElements==nzc);
1049    matrixByRow.appendRows(n_con,rowStart,column,element);
1050    delete [] rowStart;
1051    delete [] column;
1052    delete [] element;
1053    numberRows=n_con;
1054    numberColumns=n_var;
1055    //numberElements=nzc;
1056    numberBinary=nbv;
1057    numberIntegers=niv;
1058    numberAllNonLinearBoth=nlvb;
1059    numberIntegerNonLinearBoth=nlvbi;
1060    numberAllNonLinearConstraints=nlvc;
1061    numberIntegerNonLinearConstraints=nlvci;
1062    numberAllNonLinearObjective=nlvo;
1063    numberIntegerNonLinearObjective=nlvoi;
1064    /* say we want primal solution */
1065    want_xpi0=1;
1066    //double * dualSolution=NULL;
1067    // save asl
1068    // Fix memory leak one day
1069    CbcAmplInfo * info = new CbcAmplInfo;
1070    //amplGamsData_ = info;
1071    info->asl_=NULL; // as wrong form asl;
1072    info->nz_h_full_=-1; // number of nonzeros in hessian
1073    info->objval_called_with_current_x_=false;
1074    info->nerror_=0;
1075    info->obj_sign_=direction;
1076    info->conval_called_with_current_x_=false;
1077    info->non_const_x_=NULL;
1078    info->jacval_called_with_current_x_=false;
1079    info->rowStart_ = NULL;
1080    info->column_ = NULL;
1081    info->gradient_ = NULL;
1082    info->constraintValues_ = NULL;
1083  } else if (nonLinear==2) {
1084    // General nonlinear!
1085    //ASL_pfgh* asl = (ASL_pfgh*)ASL_alloc(ASL_read_pfgh);
1086    asl = ASL_alloc(ASL_read_pfgh);
1087    nl = jac0dim(stub, (long) strlen(stub));
1088    free(stub);
1089    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
1090    /* read  model*/
1091    X0 = (double*) malloc(n_var*sizeof(double));
1092    CoinZeroN(X0,n_var);
1093    // code stolen from Ipopt
1094    int retcode = pfgh_read(nl, ASL_return_read_err | ASL_findgroups);
1095
1096    switch (retcode) {
1097      case ASL_readerr_none : {}
1098      break;
1099      case ASL_readerr_nofile : {
1100        printf( "Cannot open .nl file\n");
1101        exit(-1);
1102      }
1103      break;
1104      case ASL_readerr_nonlin : {
1105        assert(false); // this better not be an error!
1106        printf( "model involves nonlinearities (ed0read)\n");
1107        exit(-1);
1108      }
1109      break;
1110      case  ASL_readerr_argerr : {
1111        printf( "user-defined function with bad args\n");
1112        exit(-1);
1113      }
1114      break;
1115      case ASL_readerr_unavail : {
1116        printf( "user-defined function not available\n");
1117        exit(-1);
1118      }
1119      break;
1120      case ASL_readerr_corrupt : {
1121        printf( "corrupt .nl file\n");
1122        exit(-1);
1123      }
1124      break;
1125      case ASL_readerr_bug : {
1126        printf( "bug in .nl reader\n");
1127        exit(-1);
1128      }
1129      break;
1130      case ASL_readerr_CLP : {
1131        printf( "ASL error message: \"solver cannot handle CLP extensions\"\n");
1132        exit(-1);
1133      }
1134      break;
1135      default: {
1136        printf( "Unknown error in stub file read. retcode = %d\n", retcode);
1137        exit(-1);
1138      }
1139      break;
1140    }
1141
1142    // see "changes" in solvers directory of ampl code...
1143    hesset(1,0,1,0,nlc);
1144
1145    assert (n_obj==1);
1146    // find the nonzero structure for the hessian
1147    // parameters to sphsetup:
1148    int coeff_obj = 1; // coefficient of the objective fn ???
1149    int mult_supplied = 1; // multipliers will be supplied
1150    int uptri = 1; // only need the upper triangular part
1151    // save asl
1152    // Fix memory leak one day
1153    CbcAmplInfo * info = new CbcAmplInfo;
1154    moreInfo_ = (void *) info;
1155    //amplGamsData_ = info;
1156    info->asl_=(ASL_pfgh *) asl;
1157    // This is not easy to get from ampl so save
1158    info->nz_h_full_= sphsetup(-1, coeff_obj, mult_supplied, uptri);
1159    info->objval_called_with_current_x_=false;
1160    info->nerror_=0;
1161    info->obj_sign_=direction;
1162    info->conval_called_with_current_x_=false;
1163    info->non_const_x_=NULL;
1164    info->jacval_called_with_current_x_=false;
1165    // Look at nonlinear
1166    if (nzc) {
1167      n_conjac[1]=nlc; // just nonlinear
1168      int * rowStart = new int [nlc+1];
1169      info->rowStart_ = rowStart;
1170      // See how many
1171      int  current_nz = 0;
1172      for (int i=0; i<nlc; i++) {
1173        for (cgrad* cg=Cgrad[i]; cg; cg = cg->next) {
1174          current_nz++;
1175        }
1176      }
1177      // setup the structure
1178      int * column = new int [current_nz];
1179      info->column_ = column;
1180      current_nz = 0;
1181      rowStart[0]=0;
1182      for (int i=0; i<nlc; i++) {
1183        for (cgrad* cg=Cgrad[i]; cg; cg = cg->next) {
1184          cg->goff=current_nz;
1185          //iRow[cg->goff] = i ;
1186          //jCol[cg->goff] = cg->varno + 1;
1187          column[cg->goff] = cg->varno ;
1188          current_nz++;
1189        }
1190        rowStart[i+1]=current_nz;
1191      }
1192      info->gradient_ = new double [nzc];
1193      info->constraintValues_ = new double [nlc];
1194    }
1195    /* objective*/
1196    objective = (double *) malloc(n_var*sizeof(double));
1197    for (i=0;i<n_var;i++)
1198      objective[i]=0.0;
1199    if (n_obj) {
1200      for (og = Ograd[0];og;og = og->next)
1201        objective[og->varno] = og->coef;
1202    }
1203    if (objtype[0])
1204      direction=-1.0;
1205    else
1206      direction=1.0;
1207    objectiveOffset_=objconst(0);
1208    /* Column bounds*/
1209    columnLower = (double *) malloc(n_var*sizeof(double));
1210    columnUpper = (double *) malloc(n_var*sizeof(double));
1211    for (i=0;i<n_var;i++) {
1212      columnLower[i]=LUv[2*i];
1213      if (columnLower[i]<= negInfinity)
1214        columnLower[i]=-COIN_DBL_MAX;
1215      columnUpper[i]=LUv[2*i+1];
1216      if (columnUpper[i]>= Infinity)
1217        columnUpper[i]=COIN_DBL_MAX;
1218    }
1219    // Build by row from scratch
1220    //matrixByRow.reserve(n_var,nzc,true);
1221    // say row orderded
1222    matrixByRow.transpose();
1223    CoinBigIndex * rowStart = new CoinBigIndex [n_con+1];
1224    int * column = new int [nzc];
1225    double * element = new double [nzc];
1226    rowStart[0]=0;
1227    numberElements=0;
1228    /* Row bounds*/
1229    rowLower = (double *) malloc(n_con*sizeof(double));
1230    rowUpper = (double *) malloc(n_con*sizeof(double));
1231    for (i=0;i<n_con;i++) {
1232      rowLower[i]=LUrhs[2*i];
1233      if (rowLower[i]<= negInfinity)
1234        rowLower[i]=-COIN_DBL_MAX;
1235      rowUpper[i]=LUrhs[2*i+1];
1236      if (rowUpper[i]>= Infinity)
1237        rowUpper[i]=COIN_DBL_MAX;
1238      for(cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
1239        column[numberElements]=cg->varno;
1240        double value = cg->coef;
1241        if (!value)
1242          value = -1.2345e-29;
1243        element[numberElements++]=value;
1244      }
1245      rowStart[i+1]=numberElements;
1246    }
1247    assert (numberElements==nzc);
1248    matrixByRow.appendRows(n_con,rowStart,column,element);
1249    delete [] rowStart;
1250    delete [] column;
1251    delete [] element;
1252    numberRows=n_con;
1253    numberColumns=n_var;
1254    numberElements=nzc;
1255    numberBinary=nbv;
1256    numberIntegers=niv;
1257    numberAllNonLinearBoth=nlvb;
1258    numberIntegerNonLinearBoth=nlvbi;
1259    numberAllNonLinearConstraints=nlvc;
1260    numberIntegerNonLinearConstraints=nlvci;
1261    numberAllNonLinearObjective=nlvo;
1262    numberIntegerNonLinearObjective=nlvoi;
1263    /* say we want primal solution */
1264    want_xpi0=1;
1265    //double * dualSolution=NULL;
1266  } else {
1267    abort();
1268  }
1269  // set problem name
1270  free(problemName_);
1271  problemName_=strdup("???");
1272
1273  // Build by row from scratch
1274  const double * element = matrixByRow.getElements();
1275  const int * column = matrixByRow.getIndices();
1276  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
1277  const int * rowLength = matrixByRow.getVectorLengths();
1278  for (i=0;i<numberRows;i++) {
1279    addRow(rowLength[i],column+rowStart[i],
1280           element+rowStart[i],rowLower[i],rowUpper[i]);
1281  }
1282  // Now do column part
1283  for (i=0;i<numberColumns;i++) {
1284    setColumnBounds(i,columnLower[i],columnUpper[i]);
1285    setColumnObjective(i,objective[i]);
1286  }
1287  for ( i=numberColumns-numberBinary-numberIntegers;
1288        i<numberColumns;i++) {
1289    setColumnIsInteger(i,true);
1290  }
1291  // and non linear
1292  for (i=numberAllNonLinearBoth-numberIntegerNonLinearBoth;
1293       i<numberAllNonLinearBoth;i++) {
1294    setColumnIsInteger(i,true);
1295  }
1296  for (i=numberAllNonLinearConstraints-numberIntegerNonLinearConstraints;
1297       i<numberAllNonLinearConstraints;i++) {
1298    setColumnIsInteger(i,true);
1299  }
1300  for (i=numberAllNonLinearObjective-numberIntegerNonLinearObjective;
1301       i<numberAllNonLinearObjective;i++) {
1302    setColumnIsInteger(i,true);
1303  }
1304  free(columnLower);
1305  free(columnUpper);
1306  free(rowLower);
1307  free(rowUpper);
1308  free(objective);
1309  // do names
1310  int iRow;
1311  for (iRow=0;iRow<numberRows_;iRow++) {
1312    char name[9];
1313    sprintf(name,"r%7.7d",iRow);
1314    setRowName(iRow,name);
1315  }
1316  int iColumn;
1317  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1318    char name[9];
1319    sprintf(name,"c%7.7d",iColumn);
1320    setColumnName(iColumn,name);
1321  }
1322  if (colqp) {
1323    // add in quadratic
1324    int nz = 1 + n_con;
1325    int nOdd=0;
1326    fint ** rowqp = colqp + nz;
1327    double ** delsqp = (double **)(rowqp + nz);
1328    for (i=0;i<=n_con;i++) {
1329      int nels = z[i];
1330      if (nels) {
1331        double * element = delsqp[i];
1332        int * start = (int *) colqp[i];
1333        int * row = (int *) rowqp[i];
1334        if (!element) {
1335          // odd row - probably not quadratic
1336          nOdd++;
1337          continue;
1338        }
1339#if 0
1340        printf("%d quadratic els\n",nels);
1341        for (int j=0;j<n_var;j++) {
1342          for (int k=start[j];k<start[j+1];k++)
1343            printf("%d %d %g\n",j,row[k],element[k]);
1344        }
1345#endif
1346        if (i) {
1347          int iRow = i-1;
1348          for (int j=0;j<n_var;j++) {
1349            for (int k=start[j];k<start[j+1];k++) {
1350              int kColumn = row[k];
1351              double value = element[k];
1352              // ampl gives twice with assumed 0.5
1353              if (kColumn<j)
1354                continue;
1355              else if (kColumn==j)
1356                value *= 0.5;
1357              const char * expr = getElementAsString(iRow,j);
1358              double constant=0.0;
1359              bool linear;
1360              if (expr&&strcmp(expr,"Numeric")) {
1361                linear=false;
1362              } else {
1363                constant = getElement(iRow,j);
1364                linear=true;
1365              }
1366              char temp[1000];
1367              char temp2[30];
1368              if (value==1.0) 
1369                sprintf(temp2,"c%7.7d",kColumn);
1370              else
1371                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1372              if (linear) {
1373                if (!constant)
1374                  strcpy(temp,temp2);
1375                else if (value>0.0) 
1376                  sprintf(temp,"%g+%s",constant,temp2);
1377                else
1378                  sprintf(temp,"%g%s",constant,temp2);
1379              } else {
1380                if (value>0.0) 
1381                  sprintf(temp,"%s+%s",expr,temp2);
1382                else
1383                  sprintf(temp,"%s%s",expr,temp2);
1384              }
1385              assert (strlen(temp)<1000);
1386              setElement(iRow,j,temp);
1387              if (amplInfo->logLevel>1)
1388                printf("el for row %d column c%7.7d is %s\n",iRow,j,temp);
1389            }
1390          }
1391        } else {
1392          // objective
1393          for (int j=0;j<n_var;j++) {
1394            for (int k=start[j];k<start[j+1];k++) {
1395              int kColumn = row[k];
1396              double value = element[k];
1397              // ampl gives twice with assumed 0.5
1398              if (kColumn<j)
1399                continue;
1400              else if (kColumn==j)
1401                value *= 0.5;
1402              const char * expr = getColumnObjectiveAsString(j);
1403              double constant=0.0;
1404              bool linear;
1405              if (expr&&strcmp(expr,"Numeric")) {
1406                linear=false;
1407              } else {
1408                constant = getColumnObjective(j);
1409                linear=true;
1410              }
1411              char temp[1000];
1412              char temp2[30];
1413              if (value==1.0) 
1414                sprintf(temp2,"c%7.7d",kColumn);
1415              else
1416                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1417              if (linear) {
1418                if (!constant)
1419                  strcpy(temp,temp2);
1420                else if (value>0.0) 
1421                  sprintf(temp,"%g+%s",constant,temp2);
1422                else
1423                  sprintf(temp,"%g%s",constant,temp2);
1424              } else {
1425                if (value>0.0) 
1426                  sprintf(temp,"%s+%s",expr,temp2);
1427                else
1428                  sprintf(temp,"%s%s",expr,temp2);
1429              }
1430              assert (strlen(temp)<1000);
1431              setObjective(j,temp);
1432              if (amplInfo->logLevel>1)
1433                printf("el for objective column c%7.7d is %s\n",j,temp);
1434            }
1435          }
1436        }
1437      }
1438    }
1439    if (nOdd) {
1440      printf("%d non-linear constraints could not be converted to quadratic\n",nOdd);
1441      exit(77);
1442    }
1443  }
1444  free(colqp);
1445  free(z);
1446  // see if any sos
1447  {
1448    char *sostype;
1449    int nsosnz, *sosbeg, *sosind, * sospri;
1450    double *sosref;
1451    int nsos;
1452    int i = ASL_suf_sos_explict_free;
1453    int copri[2], **p_sospri;
1454    copri[0] = 0;
1455    copri[1] = 0;
1456    p_sospri = &sospri;
1457    nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1458                   &sosbeg, &sosind, &sosref);
1459    if (nsos) {
1460      numberSOS_=nsos;
1461      typeSOS_ = new int [numberSOS_];
1462      prioritySOS_ = new int [numberSOS_];
1463      startSOS_ = new int [numberSOS_+1];
1464      memberSOS_ = new int[nsosnz];
1465      referenceSOS_ = new double [nsosnz];
1466      sos_kludge(nsos, sosbeg, sosref,sosind);
1467      for (int i=0;i<nsos;i++) {
1468        int ichar = sostype[i];
1469        assert (ichar=='1'||ichar=='2');
1470        typeSOS_[i]=ichar-'0';
1471      } 
1472      memcpy(prioritySOS_,sospri,nsos*sizeof(int));
1473      memcpy(startSOS_,sosbeg,(nsos+1)*sizeof(int));
1474      memcpy(memberSOS_,sosind,nsosnz*sizeof(int));
1475      memcpy(referenceSOS_,sosref,nsosnz*sizeof(double));
1476    }
1477  }
1478}
1479#else
1480#include "Cbc_ampl.h"
1481int
1482readAmpl(ampl_info * info, int argc, char **argv, void ** coinModel)
1483{
1484  return 0;
1485}
1486void freeArrays1(ampl_info * info)
1487{
1488}
1489void freeArrays2(ampl_info * info)
1490{
1491}
1492void freeArgs(ampl_info * info)
1493{
1494}
1495int ampl_obj_prec()
1496{
1497  return 0;
1498}
1499void writeAmpl(ampl_info * info)
1500{
1501}
1502#endif
Note: See TracBrowser for help on using the repository browser.