source: stable/2.4/Cbc/src/Cbc_ampl.cpp @ 1271

Last change on this file since 1271 was 1271, checked in by forrest, 10 years ago

Creating new stable branch 2.4 from trunk (rev 1270)

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