source: stable/2.3/Cbc/src/Cbc_ampl.cpp @ 1142

Last change on this file since 1142 was 1125, checked in by forrest, 11 years ago

ampl compiler errors

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