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

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

Added extra return at end of each source file where needed, to remove possible linefeed conflicts (NightlyBuild? errors)

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