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

Last change on this file since 1733 was 1643, checked in by stefan, 8 years ago

apply fixes for compiler warnings send by Alpar Juettner:

  • Every appearance of DBL_MAX was replaced by COIN_DBL_MAX.
  • The unused parameter warnings were resolved by removing their names from the declaration.
  • Finally, I weeded out all the old-style casts from the code.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 52.8 KB
Line 
1/* $Id: Cbc_ampl.cpp 1643 2011-04-29 18:06:50Z stefan $ */
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        for (i = 0; i < n_var; i++) {
474            columnLower[i] = LUv[2*i];
475            if (columnLower[i] <= negInfinity)
476                columnLower[i] = -COIN_DBL_MAX;
477            columnUpper[i] = LUv[2*i+1];
478            if (columnUpper[i] >= Infinity)
479                columnUpper[i] = COIN_DBL_MAX;
480        }
481        /* Row bounds*/
482        rowLower = (double *) malloc(n_con * sizeof(double));
483        rowUpper = (double *) malloc(n_con * sizeof(double));
484        for (i = 0; i < n_con; i++) {
485            rowLower[i] = LUrhs[2*i];
486            if (rowLower[i] <= negInfinity)
487                rowLower[i] = -COIN_DBL_MAX;
488            rowUpper[i] = LUrhs[2*i+1];
489            if (rowUpper[i] >= Infinity)
490                rowUpper[i] = COIN_DBL_MAX;
491        }
492        info->numberRows = n_con;
493        info->numberColumns = n_var;
494        info->numberElements = nzc;
495        info->numberBinary = nbv;
496        info->numberIntegers = niv + nbv;
497        info->objective = obj;
498        info->rowLower = rowLower;
499        info->rowUpper = rowUpper;
500        info->columnLower = columnLower;
501        info->columnUpper = columnUpper;
502        info->starts = A_colstarts;
503        /*A_colstarts=NULL;*/
504        info->rows = A_rownos;
505        /*A_rownos=NULL;*/
506        info->elements = A_vals;
507        /*A_vals=NULL;*/
508        info->primalSolution = NULL;
509        /* put in primalSolution if exists */
510        if (X0) {
511            info->primalSolution = (double *) malloc(n_var * sizeof(double));
512            memcpy(info->primalSolution, X0, n_var*sizeof(double));
513        }
514        info->dualSolution = NULL;
515        if (niv + nbv > 0)
516            mip_stuff(); // get any extra info
517        if ((!(niv + nbv) && (csd->kind & ASL_Sufkind_input))
518                || (rsd->kind & ASL_Sufkind_input)) {
519            /* convert status - need info on map */
520            static int map[] = {1, 3, 1, 1, 2, 1, 1};
521            stat_map(info->columnStatus, n_var, map, 6, "incoming columnStatus");
522            stat_map(info->rowStatus, n_con, map, 6, "incoming rowStatus");
523        } else {
524            /* all slack basis */
525            // leave status for output */
526#ifdef JJF_ZERO
527            free(info->rowStatus);
528            info->rowStatus = NULL;
529            free(info->columnStatus);
530            info->columnStatus = NULL;
531#endif
532        }
533    } else {
534        // QP
535        // Add .nl if not there
536        if (!strstr(fileName, ".nl"))
537            strcat(fileName, ".nl");
538        CoinModel * model = new CoinModel((nonLinearType > 10) ? 2 : 1, fileName, info);
539        if (model->numberRows() > 0 || model->numberColumns() > 0)
540            *coinModel = (void *) model;
541        Oinfo.uinfo = tempBuffer;
542        if (getopts(argv, &Oinfo))
543            return 1;
544        Oinfo.wantsol = 1;
545        if (objtype[0])
546            info->direction = -1.0;
547        else
548            info->direction = 1.0;
549        model->setOptimizationDirection(info->direction);
550        info->offset = objconst(0);
551        info->numberRows = n_con;
552        info->numberColumns = n_var;
553        info->numberElements = nzc;
554        info->numberBinary = nbv;
555        int numberIntegers = niv + nlvci + nlvoi + nbv;
556        if (nlvci + nlvoi + nlvc + nlvo) {
557            // Non linear
558            // No idea if there are overlaps so compute
559            int numberIntegers = 0;
560            for ( i = 0; i < n_var; i++) {
561                if (model->columnIsInteger(i))
562                    numberIntegers++;
563            }
564        }
565        info->numberIntegers = numberIntegers;
566        // Say nonlinear if it is
567        info->nonLinear = nlvc + nlvo;
568        if (numberIntegers > 0) {
569            mip_stuff(); // get any extra info
570            if (info->cut)
571                model->setCutMarker(info->numberRows, info->cut);
572            if (info->priorities)
573                model->setPriorities(info->numberColumns, info->priorities);
574        }
575    }
576    /* add -solve - unless something there already
577     - also check for sleep=yes */
578    {
579        int found = 0;
580        int foundLog = 0;
581        int foundSleep = 0;
582        const char * something[] = {"solve", "branch", "duals", "primals", "user"};
583        for (i = 0; i < info->numberArguments; i++) {
584            unsigned int j;
585            const char * argument = info->arguments[i];
586            for (j = 0; j < sizeof(something) / sizeof(char *); j++) {
587                const char * check = something[j];
588                if (!strncmp(argument, check, sizeof(check))) {
589                    found = (int)(j + 1);
590                } else if (!strncmp(argument, "log", 3)) {
591                    foundLog = 1;
592                } else if (!strncmp(argument, "sleep", 5)) {
593                    foundSleep = 1;
594                }
595            }
596        }
597        if (foundLog) {
598            /* print options etc */
599            for (i = 0; i < saveArgc; i++)
600                printf("%s ", saveArgv[i]);
601            printf("\n");
602            if (environment)
603                printf("env %s\n", environment);
604            /*printf("%d rows %d columns %d elements\n",n_con,n_var,nzc);*/
605        }
606        if (!found) {
607            if (!strlen(algFound)) {
608                info->arguments = (char **) realloc(info->arguments, (info->numberArguments + 1) * sizeof(char *));
609                info->arguments[info->numberArguments++] = strdup("-solve");
610            } else {
611                // use algorithm from keyword
612                info->arguments = (char **) realloc(info->arguments, (info->numberArguments + 1) * sizeof(char *));
613                info->arguments[info->numberArguments++] = strdup(algFound);
614            }
615        }
616        if (foundSleep) {
617            /* let user copy .nl file */
618            fprintf(stderr, "You can copy .nl file %s for debug purposes or attach debugger\n", saveArgv[1]);
619            fprintf(stderr, "Type q to quit, anything else to continue\n");
620            char getChar = getc(stdin);
621            if (getChar == 'q' || getChar == 'Q')
622                exit(1);
623        }
624    }
625    /* add -quit */
626    info->arguments = (char **) realloc(info->arguments, (info->numberArguments + 1) * sizeof(char *));
627    info->arguments[info->numberArguments++] = strdup("-quit");
628    return 0;
629}
630void freeArrays1(ampl_info * info)
631{
632    free(info->objective);
633    info->objective = NULL;
634    free(info->rowLower);
635    info->rowLower = NULL;
636    free(info->rowUpper);
637    info->rowUpper = NULL;
638    free(info->columnLower);
639    info->columnLower = NULL;
640    free(info->columnUpper);
641    info->columnUpper = NULL;
642    /* this one not freed by ASL_free */
643    free(info->elements);
644    info->elements = NULL;
645    free(info->primalSolution);
646    info->primalSolution = NULL;
647    free(info->dualSolution);
648    info->dualSolution = NULL;
649    /*free(info->rowStatus);
650    info->rowStatus=NULL;
651    free(info->columnStatus);
652    info->columnStatus=NULL;*/
653}
654void freeArrays2(ampl_info * info)
655{
656    free(info->primalSolution);
657    info->primalSolution = NULL;
658    free(info->dualSolution);
659    info->dualSolution = NULL;
660    free(info->rowStatus);
661    info->rowStatus = NULL;
662    free(info->columnStatus);
663    info->columnStatus = NULL;
664    free(info->priorities);
665    info->priorities = NULL;
666    free(info->branchDirection);
667    info->branchDirection = NULL;
668    free(info->pseudoDown);
669    info->pseudoDown = NULL;
670    free(info->pseudoUp);
671    info->pseudoUp = NULL;
672    free(info->sosType);
673    info->sosType = NULL;
674    free(info->sosPriority);
675    info->sosPriority = NULL;
676    free(info->sosStart);
677    info->sosStart = NULL;
678    free(info->sosIndices);
679    info->sosIndices = NULL;
680    free(info->sosReference);
681    info->sosReference = NULL;
682    free(info->cut);
683    info->cut = NULL;
684    ASL_free(&asl);
685}
686void freeArgs(ampl_info * info)
687{
688    int i;
689    for ( i = 0; i < info->numberArguments; i++)
690        free(info->arguments[i]);
691    free(info->arguments);
692}
693int ampl_obj_prec()
694{
695    return obj_prec();
696}
697void writeAmpl(ampl_info * info)
698{
699    char buf[1000];
700    typedef struct {
701        const char *msg;
702        int code;
703        int wantObj;
704    } Sol_info;
705    static Sol_info solinfo[] = {
706        { "optimal solution",                   000, 1 },
707        { "infeasible",                         200, 1 },
708        { "unbounded",                          300, 0 },
709        { "iteration limit etc",                        400, 1 },
710        { "solution limit",                             401, 1 },
711        { "ran out of space",                   500, 0 },
712        { "status unknown",                             501, 1 },
713        { "bug!",                                       502, 0 },
714        { "best MIP solution so far restored",  101, 1 },
715        { "failed to restore best MIP solution",        503, 1 },
716        { "optimal (?) solution",                       100, 1 }
717    };
718    /* convert status - need info on map */
719    static int map[] = {0, 3, 4, 1};
720    sprintf(buf, "%s %s", Oinfo.bsname, info->buffer);
721    solve_result_num = solinfo[info->problemStatus].code;
722    if (info->columnStatus) {
723        stat_map(info->columnStatus, n_var, map, 4, "outgoing columnStatus");
724        stat_map(info->rowStatus, n_con, map, 4, "outgoing rowStatus");
725        suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
726        suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
727    }
728    write_sol(buf, info->primalSolution, info->dualSolution, &Oinfo);
729}
730/* Read a problem from AMPL nl file
731 */
732CoinModel::CoinModel( int nonLinear, const char * fileName, const void * info)
733        :  CoinBaseModel(),
734        maximumRows_(0),
735        maximumColumns_(0),
736        numberElements_(0),
737        maximumElements_(0),
738        numberQuadraticElements_(0),
739        maximumQuadraticElements_(0),
740        rowLower_(NULL),
741        rowUpper_(NULL),
742        rowType_(NULL),
743        objective_(NULL),
744        columnLower_(NULL),
745        columnUpper_(NULL),
746        integerType_(NULL),
747        columnType_(NULL),
748        start_(NULL),
749        elements_(NULL),
750        packedMatrix_(NULL),
751        quadraticElements_(NULL),
752        sortIndices_(NULL),
753        sortElements_(NULL),
754        sortSize_(0),
755        sizeAssociated_(0),
756        associated_(NULL),
757        numberSOS_(0),
758        startSOS_(NULL),
759        memberSOS_(NULL),
760        typeSOS_(NULL),
761        prioritySOS_(NULL),
762        referenceSOS_(NULL),
763        priority_(NULL),
764        cut_(NULL),
765        moreInfo_(NULL),
766        type_(-1),
767        links_(0)
768{
769    problemName_ = "";
770    int status = 0;
771    if (!strcmp(fileName, "-") || !strcmp(fileName, "stdin")) {
772        // stdin
773    } else {
774        std::string name = fileName;
775        bool readable = fileCoinReadable(name);
776        if (!readable) {
777            std::cerr << "Unable to open file "
778                      << fileName << std::endl;
779            status = -1;
780        }
781    }
782    if (!status) {
783        gdb(nonLinear, fileName, info);
784    }
785}
786#ifdef JJF_ZERO
787static real
788qterm(ASL *asl, fint *colq, fint *rowq, real *delsq)
789{
790    double t, t1, *x, *x0, *xe;
791    fint *rq0, *rqe;
792
793    t = 0.;
794    x0 = x = X0;
795    xe = x + n_var;
796    rq0 = rowq;
797    while (x < xe) {
798        t1 = *x++;
799        rqe = rq0 + *++colq;
800        while (rowq < rqe)
801            t += t1 * x0[*rowq++]**delsq++;
802    }
803    return 0.5 * t;
804}
805#endif
806// stolen from IPopt with changes
807typedef struct {
808    double obj_sign_;
809    ASL_pfgh * asl_;
810    double * non_const_x_;
811    int * column_; // for jacobian
812    int * rowStart_;
813    double * gradient_;
814    double * constraintValues_;
815    int nz_h_full_; // number of nonzeros in hessian
816    int nerror_;
817    bool objval_called_with_current_x_;
818    bool conval_called_with_current_x_;
819    bool jacval_called_with_current_x_;
820} CbcAmplInfo;
821
822void
823CoinModel::gdb( int nonLinear, const char * fileName, const void * info)
824{
825    const ampl_info * amplInfo = (const ampl_info *) info;
826    ograd *og = NULL;
827    int i;
828    SufDesc *csd = NULL;
829    SufDesc *rsd = NULL;
830    /*bool *basis, *lower;*/
831    /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
832    //char tempBuffer[20];
833    double * objective = NULL;
834    double * columnLower = NULL;
835    double * columnUpper = NULL;
836    double * rowLower = NULL;
837    double * rowUpper = NULL;
838    int * columnStatus = NULL;
839    int * rowStatus = NULL;
840    int numberRows = -1;
841    int numberColumns = -1;
842    int numberElements = -1;
843    int numberBinary = -1;
844    int numberIntegers = -1;
845    int numberAllNonLinearBoth = 0;
846    int numberIntegerNonLinearBoth = 0;
847    int numberAllNonLinearConstraints = 0;
848    int numberIntegerNonLinearConstraints = 0;
849    int numberAllNonLinearObjective = 0;
850    int numberIntegerNonLinearObjective = 0;
851    double * primalSolution = NULL;
852    double direction = 1.0;
853    char * stub = strdup(fileName);
854    CoinPackedMatrix matrixByRow;
855    fint ** colqp = NULL;
856    int *z = NULL;
857    if (nonLinear == 0) {
858        // linear
859        asl = ASL_alloc(ASL_read_f);
860        nl = jac0dim(stub, 0);
861        free(stub);
862        suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
863
864        /* set A_vals to get the constraints column-wise (malloc so can be freed) */
865        A_vals = (double *) malloc(nzc * sizeof(double));
866        if (!A_vals) {
867            printf("no memory\n");
868            return ;
869        }
870        /* say we want primal solution */
871        want_xpi0 = 1;
872        /* for basis info */
873        columnStatus = (int *) malloc(n_var * sizeof(int));
874        rowStatus = (int *) malloc(n_con * sizeof(int));
875        csd = suf_iput("sstatus", ASL_Sufkind_var, columnStatus);
876        rsd = suf_iput("sstatus", ASL_Sufkind_con, rowStatus);
877        /* read linear model*/
878        f_read(nl, 0);
879        // see if any sos
880        if (true) {
881            char *sostype;
882            int nsosnz, *sosbeg, *sosind, * sospri;
883            double *sosref;
884            int nsos;
885            int i = ASL_suf_sos_explict_free;
886            int copri[2], **p_sospri;
887            copri[0] = 0;
888            copri[1] = 0;
889            p_sospri = &sospri;
890            nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
891                           &sosbeg, &sosind, &sosref);
892            if (nsos) {
893                abort();
894#ifdef JJF_ZERO
895                info->numberSos = nsos;
896                info->sosType = (char *) malloc(nsos);
897                info->sosPriority = (int *) malloc(nsos * sizeof(int));
898                info->sosStart = (int *) malloc((nsos + 1) * sizeof(int));
899                info->sosIndices = (int *) malloc(nsosnz * sizeof(int));
900                info->sosReference = (double *) malloc(nsosnz * sizeof(double));
901                sos_kludge(nsos, sosbeg, sosref, sosind);
902                for (int i = 0; i < nsos; i++) {
903                    int ichar = sostype[i];
904                    assert (ichar == '1' || ichar == '2');
905                    info->sosType[i] = ichar - '0';
906                }
907                memcpy(info->sosPriority, sospri, nsos*sizeof(int));
908                memcpy(info->sosStart, sosbeg, (nsos + 1)*sizeof(int));
909                memcpy(info->sosIndices, sosind, nsosnz*sizeof(int));
910                memcpy(info->sosReference, sosref, nsosnz*sizeof(double));
911#endif
912            }
913        }
914
915        /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
916        //Oinfo.uinfo = tempBuffer;
917        //if (getopts(argv, &Oinfo))
918        //return 1;
919        /* objective*/
920        objective = (double *) malloc(n_var * sizeof(double));
921        for (i = 0; i < n_var; i++)
922            objective[i] = 0.0;
923        if (n_obj) {
924            for (og = Ograd[0]; og; og = og->next)
925                objective[og->varno] = og->coef;
926        }
927        if (objtype[0])
928            direction = -1.0;
929        else
930            direction = 1.0;
931        objectiveOffset_ = objconst(0);
932        /* Column bounds*/
933        columnLower = (double *) malloc(n_var * sizeof(double));
934        columnUpper = (double *) malloc(n_var * sizeof(double));
935        for (i = 0; i < n_var; i++) {
936            columnLower[i] = LUv[2*i];
937            if (columnLower[i] <= negInfinity)
938                columnLower[i] = -COIN_DBL_MAX;
939            columnUpper[i] = LUv[2*i+1];
940            if (columnUpper[i] >= Infinity)
941                columnUpper[i] = COIN_DBL_MAX;
942        }
943        /* Row bounds*/
944        rowLower = (double *) malloc(n_con * sizeof(double));
945        rowUpper = (double *) malloc(n_con * sizeof(double));
946        for (i = 0; i < n_con; i++) {
947            rowLower[i] = LUrhs[2*i];
948            if (rowLower[i] <= negInfinity)
949                rowLower[i] = -COIN_DBL_MAX;
950            rowUpper[i] = LUrhs[2*i+1];
951            if (rowUpper[i] >= Infinity)
952                rowUpper[i] = COIN_DBL_MAX;
953        }
954        numberRows = n_con;
955        numberColumns = n_var;
956        numberElements = nzc;
957        numberBinary = nbv;
958        numberIntegers = niv;
959        /* put in primalSolution if exists */
960        if (X0) {
961            primalSolution = (double *) malloc(n_var * sizeof(double));
962            memcpy( primalSolution, X0, n_var*sizeof(double));
963        }
964        //double * dualSolution=NULL;
965        if (niv + nbv > 0)
966            mip_stuff(); // get any extra info
967        if ((!(niv + nbv) && (csd->kind & ASL_Sufkind_input))
968                || (rsd->kind & ASL_Sufkind_input)) {
969            /* convert status - need info on map */
970            static int map[] = {1, 3, 1, 1, 2, 1, 1};
971            stat_map(columnStatus, n_var, map, 6, "incoming columnStatus");
972            stat_map(rowStatus, n_con, map, 6, "incoming rowStatus");
973        } else {
974            /* all slack basis */
975            // leave status for output */
976#ifdef JJF_ZERO
977            free(rowStatus);
978            rowStatus = NULL;
979            free(columnStatus);
980            columnStatus = NULL;
981#endif
982        }
983        CoinPackedMatrix columnCopy(true, numberRows, numberColumns, numberElements,
984                                    A_vals, A_rownos, A_colstarts, NULL);
985        matrixByRow.reverseOrderedCopyOf(columnCopy);
986    } else if (nonLinear == 1) {
987        // quadratic
988        asl = ASL_alloc(ASL_read_fg);
989        nl = jac0dim(stub, (long) strlen(stub));
990        free(stub);
991        suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
992        /* read  model*/
993        X0 = (double*) malloc(n_var * sizeof(double));
994        CoinZeroN(X0, n_var);
995        qp_read(nl, 0);
996        assert (n_obj == 1);
997        int nz = 1 + n_con;
998        colqp = (fint**) malloc(nz * (2 * sizeof(int*)
999                                      + sizeof(double*)));
1000        fint ** rowqp = colqp + nz;
1001        double ** delsqp = (double **)(rowqp + nz);
1002        z = (int*) malloc(nz * sizeof(int));
1003        for (i = 0; i <= n_con; i++) {
1004            z[i] = nqpcheck(-i, rowqp + i, colqp + i, delsqp + i);
1005        }
1006        qp_opify();
1007        /* objective*/
1008        objective = (double *) malloc(n_var * sizeof(double));
1009        for (i = 0; i < n_var; i++)
1010            objective[i] = 0.0;
1011        if (n_obj) {
1012            for (og = Ograd[0]; og; og = og->next)
1013                objective[og->varno] = og->coef;
1014        }
1015        if (objtype[0])
1016            direction = -1.0;
1017        else
1018            direction = 1.0;
1019        objectiveOffset_ = objconst(0);
1020        /* Column bounds*/
1021        columnLower = (double *) malloc(n_var * sizeof(double));
1022        columnUpper = (double *) malloc(n_var * sizeof(double));
1023        for (i = 0; i < n_var; i++) {
1024            columnLower[i] = LUv[2*i];
1025            if (columnLower[i] <= negInfinity)
1026                columnLower[i] = -COIN_DBL_MAX;
1027            columnUpper[i] = LUv[2*i+1];
1028            if (columnUpper[i] >= Infinity)
1029                columnUpper[i] = COIN_DBL_MAX;
1030        }
1031        // Build by row from scratch
1032        //matrixByRow.reserve(n_var,nzc,true);
1033        // say row orderded
1034        matrixByRow.transpose();
1035        /* Row bounds*/
1036        rowLower = (double *) malloc(n_con * sizeof(double));
1037        rowUpper = (double *) malloc(n_con * sizeof(double));
1038        CoinBigIndex * rowStart = new CoinBigIndex [n_con+1];
1039        int * column = new int [nzc];
1040        double * element = new double [nzc];
1041        rowStart[0] = 0;
1042        numberElements = 0;
1043        for (i = 0; i < n_con; i++) {
1044            rowLower[i] = LUrhs[2*i];
1045            if (rowLower[i] <= negInfinity)
1046                rowLower[i] = -COIN_DBL_MAX;
1047            rowUpper[i] = LUrhs[2*i+1];
1048            if (rowUpper[i] >= Infinity)
1049                rowUpper[i] = COIN_DBL_MAX;
1050            for (cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
1051                column[numberElements] = cg->varno;
1052                element[numberElements++] = cg->coef;
1053            }
1054            rowStart[i+1] = numberElements;
1055        }
1056        assert (numberElements == nzc);
1057        matrixByRow.appendRows(n_con, rowStart, column, element);
1058        delete [] rowStart;
1059        delete [] column;
1060        delete [] element;
1061        numberRows = n_con;
1062        numberColumns = n_var;
1063        //numberElements=nzc;
1064        numberBinary = nbv;
1065        numberIntegers = niv;
1066        numberAllNonLinearBoth = nlvb;
1067        numberIntegerNonLinearBoth = nlvbi;
1068        numberAllNonLinearConstraints = nlvc;
1069        numberIntegerNonLinearConstraints = nlvci;
1070        numberAllNonLinearObjective = nlvo;
1071        numberIntegerNonLinearObjective = nlvoi;
1072        /* say we want primal solution */
1073        want_xpi0 = 1;
1074        //double * dualSolution=NULL;
1075        // save asl
1076        // Fix memory leak one day
1077        CbcAmplInfo * info = new CbcAmplInfo;
1078        //amplGamsData_ = info;
1079        info->asl_ = NULL; // as wrong form asl;
1080        info->nz_h_full_ = -1; // number of nonzeros in hessian
1081        info->objval_called_with_current_x_ = false;
1082        info->nerror_ = 0;
1083        info->obj_sign_ = direction;
1084        info->conval_called_with_current_x_ = false;
1085        info->non_const_x_ = NULL;
1086        info->jacval_called_with_current_x_ = false;
1087        info->rowStart_ = NULL;
1088        info->column_ = NULL;
1089        info->gradient_ = NULL;
1090        info->constraintValues_ = NULL;
1091    } else if (nonLinear == 2) {
1092        // General nonlinear!
1093        //ASL_pfgh* asl = (ASL_pfgh*)ASL_alloc(ASL_read_pfgh);
1094        asl = ASL_alloc(ASL_read_pfgh);
1095        nl = jac0dim(stub, (long) strlen(stub));
1096        free(stub);
1097        suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
1098        /* read  model*/
1099        X0 = (double*) malloc(n_var * sizeof(double));
1100        CoinZeroN(X0, n_var);
1101        // code stolen from Ipopt
1102        int retcode = pfgh_read(nl, ASL_return_read_err | ASL_findgroups);
1103
1104        switch (retcode) {
1105        case ASL_readerr_none : {}
1106        break;
1107        case ASL_readerr_nofile : {
1108            printf( "Cannot open .nl file\n");
1109            exit(-1);
1110        }
1111        break;
1112        case ASL_readerr_nonlin : {
1113            assert(false); // this better not be an error!
1114            printf( "model involves nonlinearities (ed0read)\n");
1115            exit(-1);
1116        }
1117        break;
1118        case  ASL_readerr_argerr : {
1119            printf( "user-defined function with bad args\n");
1120            exit(-1);
1121        }
1122        break;
1123        case ASL_readerr_unavail : {
1124            printf( "user-defined function not available\n");
1125            exit(-1);
1126        }
1127        break;
1128        case ASL_readerr_corrupt : {
1129            printf( "corrupt .nl file\n");
1130            exit(-1);
1131        }
1132        break;
1133        case ASL_readerr_bug : {
1134            printf( "bug in .nl reader\n");
1135            exit(-1);
1136        }
1137        break;
1138        case ASL_readerr_CLP : {
1139            printf( "ASL error message: \"solver cannot handle CLP extensions\"\n");
1140            exit(-1);
1141        }
1142        break;
1143        default: {
1144            printf( "Unknown error in stub file read. retcode = %d\n", retcode);
1145            exit(-1);
1146        }
1147        break;
1148        }
1149
1150        // see "changes" in solvers directory of ampl code...
1151        hesset(1, 0, 1, 0, nlc);
1152
1153        assert (n_obj == 1);
1154        // find the nonzero structure for the hessian
1155        // parameters to sphsetup:
1156        int coeff_obj = 1; // coefficient of the objective fn ???
1157        int mult_supplied = 1; // multipliers will be supplied
1158        int uptri = 1; // only need the upper triangular part
1159        // save asl
1160        // Fix memory leak one day
1161        CbcAmplInfo * info = new CbcAmplInfo;
1162        moreInfo_ = (void *) info;
1163        //amplGamsData_ = info;
1164        info->asl_ = (ASL_pfgh *) asl;
1165        // This is not easy to get from ampl so save
1166        info->nz_h_full_ = sphsetup(-1, coeff_obj, mult_supplied, uptri);
1167        info->objval_called_with_current_x_ = false;
1168        info->nerror_ = 0;
1169        info->obj_sign_ = direction;
1170        info->conval_called_with_current_x_ = false;
1171        info->non_const_x_ = NULL;
1172        info->jacval_called_with_current_x_ = false;
1173        // Look at nonlinear
1174        if (nzc) {
1175            n_conjac[1] = nlc; // just nonlinear
1176            int * rowStart = new int [nlc+1];
1177            info->rowStart_ = rowStart;
1178            // See how many
1179            int  current_nz = 0;
1180            for (int i = 0; i < nlc; i++) {
1181                for (cgrad* cg = Cgrad[i]; cg; cg = cg->next) {
1182                    current_nz++;
1183                }
1184            }
1185            // setup the structure
1186            int * column = new int [current_nz];
1187            info->column_ = column;
1188            current_nz = 0;
1189            rowStart[0] = 0;
1190            for (int i = 0; i < nlc; i++) {
1191                for (cgrad* cg = Cgrad[i]; cg; cg = cg->next) {
1192                    cg->goff = current_nz;
1193                    //iRow[cg->goff] = i ;
1194                    //jCol[cg->goff] = cg->varno + 1;
1195                    column[cg->goff] = cg->varno ;
1196                    current_nz++;
1197                }
1198                rowStart[i+1] = current_nz;
1199            }
1200            info->gradient_ = new double [nzc];
1201            info->constraintValues_ = new double [nlc];
1202        }
1203        /* objective*/
1204        objective = (double *) malloc(n_var * sizeof(double));
1205        for (i = 0; i < n_var; i++)
1206            objective[i] = 0.0;
1207        if (n_obj) {
1208            for (og = Ograd[0]; og; og = og->next)
1209                objective[og->varno] = og->coef;
1210        }
1211        if (objtype[0])
1212            direction = -1.0;
1213        else
1214            direction = 1.0;
1215        objectiveOffset_ = objconst(0);
1216        /* Column bounds*/
1217        columnLower = (double *) malloc(n_var * sizeof(double));
1218        columnUpper = (double *) malloc(n_var * sizeof(double));
1219        for (i = 0; i < n_var; i++) {
1220            columnLower[i] = LUv[2*i];
1221            if (columnLower[i] <= negInfinity)
1222                columnLower[i] = -COIN_DBL_MAX;
1223            columnUpper[i] = LUv[2*i+1];
1224            if (columnUpper[i] >= Infinity)
1225                columnUpper[i] = COIN_DBL_MAX;
1226        }
1227        // Build by row from scratch
1228        //matrixByRow.reserve(n_var,nzc,true);
1229        // say row orderded
1230        matrixByRow.transpose();
1231        CoinBigIndex * rowStart = new CoinBigIndex [n_con+1];
1232        int * column = new int [nzc];
1233        double * element = new double [nzc];
1234        rowStart[0] = 0;
1235        numberElements = 0;
1236        /* Row bounds*/
1237        rowLower = (double *) malloc(n_con * sizeof(double));
1238        rowUpper = (double *) malloc(n_con * sizeof(double));
1239        for (i = 0; i < n_con; i++) {
1240            rowLower[i] = LUrhs[2*i];
1241            if (rowLower[i] <= negInfinity)
1242                rowLower[i] = -COIN_DBL_MAX;
1243            rowUpper[i] = LUrhs[2*i+1];
1244            if (rowUpper[i] >= Infinity)
1245                rowUpper[i] = COIN_DBL_MAX;
1246            for (cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
1247                column[numberElements] = cg->varno;
1248                double value = cg->coef;
1249                if (!value)
1250                    value = -1.2345e-29;
1251                element[numberElements++] = value;
1252            }
1253            rowStart[i+1] = numberElements;
1254        }
1255        assert (numberElements == nzc);
1256        matrixByRow.appendRows(n_con, rowStart, column, element);
1257        delete [] rowStart;
1258        delete [] column;
1259        delete [] element;
1260        numberRows = n_con;
1261        numberColumns = n_var;
1262        numberElements = nzc;
1263        numberBinary = nbv;
1264        numberIntegers = niv;
1265        numberAllNonLinearBoth = nlvb;
1266        numberIntegerNonLinearBoth = nlvbi;
1267        numberAllNonLinearConstraints = nlvc;
1268        numberIntegerNonLinearConstraints = nlvci;
1269        numberAllNonLinearObjective = nlvo;
1270        numberIntegerNonLinearObjective = nlvoi;
1271        /* say we want primal solution */
1272        want_xpi0 = 1;
1273        //double * dualSolution=NULL;
1274    } else {
1275        abort();
1276    }
1277    // set problem name
1278    problemName_ = "???";
1279
1280    // Build by row from scratch
1281    const double * element = matrixByRow.getElements();
1282    const int * column = matrixByRow.getIndices();
1283    const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
1284    const int * rowLength = matrixByRow.getVectorLengths();
1285    for (i = 0; i < numberRows; i++) {
1286        addRow(rowLength[i], column + rowStart[i],
1287               element + rowStart[i], rowLower[i], rowUpper[i]);
1288    }
1289    // Now do column part
1290    for (i = 0; i < numberColumns; i++) {
1291        setColumnBounds(i, columnLower[i], columnUpper[i]);
1292        setColumnObjective(i, objective[i]);
1293    }
1294    for ( i = numberColumns - numberBinary - numberIntegers;
1295            i < numberColumns; i++) {
1296        setColumnIsInteger(i, true);
1297    }
1298    // and non linear
1299    for (i = numberAllNonLinearBoth - numberIntegerNonLinearBoth;
1300            i < numberAllNonLinearBoth; i++) {
1301        setColumnIsInteger(i, true);
1302    }
1303    for (i = numberAllNonLinearConstraints - numberIntegerNonLinearConstraints;
1304            i < numberAllNonLinearConstraints; i++) {
1305        setColumnIsInteger(i, true);
1306    }
1307    for (i = numberAllNonLinearObjective - numberIntegerNonLinearObjective;
1308            i < numberAllNonLinearObjective; i++) {
1309        setColumnIsInteger(i, true);
1310    }
1311    free(columnLower);
1312    free(columnUpper);
1313    free(rowLower);
1314    free(rowUpper);
1315    free(objective);
1316    // do names
1317    int iRow;
1318    for (iRow = 0; iRow < numberRows_; iRow++) {
1319        char name[9];
1320        sprintf(name, "r%7.7d", iRow);
1321        setRowName(iRow, name);
1322    }
1323    int iColumn;
1324    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1325        char name[9];
1326        sprintf(name, "c%7.7d", iColumn);
1327        setColumnName(iColumn, name);
1328    }
1329    if (colqp) {
1330        // add in quadratic
1331        int nz = 1 + n_con;
1332        int nOdd = 0;
1333        fint ** rowqp = colqp + nz;
1334        double ** delsqp = (double **)(rowqp + nz);
1335        for (i = 0; i <= n_con; i++) {
1336            int nels = z[i];
1337            if (nels) {
1338                double * element = delsqp[i];
1339                int * start = (int *) colqp[i];
1340                int * row = (int *) rowqp[i];
1341                if (!element) {
1342                    // odd row - probably not quadratic
1343                    nOdd++;
1344                    continue;
1345                }
1346#ifdef JJF_ZERO
1347                printf("%d quadratic els\n", nels);
1348                for (int j = 0; j < n_var; j++) {
1349                    for (int k = start[j]; k < start[j+1]; k++)
1350                        printf("%d %d %g\n", j, row[k], element[k]);
1351                }
1352#endif
1353                if (i) {
1354                    int iRow = i - 1;
1355                    for (int j = 0; j < n_var; j++) {
1356                        for (int k = start[j]; k < start[j+1]; k++) {
1357                            int kColumn = row[k];
1358                            double value = element[k];
1359                            // ampl gives twice with assumed 0.5
1360                            if (kColumn < j)
1361                                continue;
1362                            else if (kColumn == j)
1363                                value *= 0.5;
1364                            const char * expr = getElementAsString(iRow, j);
1365                            double constant = 0.0;
1366                            bool linear;
1367                            if (expr && strcmp(expr, "Numeric")) {
1368                                linear = false;
1369                            } else {
1370                                constant = getElement(iRow, j);
1371                                linear = true;
1372                            }
1373                            char temp[1000];
1374                            char temp2[30];
1375                            if (value == 1.0)
1376                                sprintf(temp2, "c%7.7d", kColumn);
1377                            else
1378                                sprintf(temp2, "%g*c%7.7d", value, kColumn);
1379                            if (linear) {
1380                                if (!constant)
1381                                    strcpy(temp, temp2);
1382                                else if (value > 0.0)
1383                                    sprintf(temp, "%g+%s", constant, temp2);
1384                                else
1385                                    sprintf(temp, "%g%s", constant, temp2);
1386                            } else {
1387                                if (value > 0.0)
1388                                    sprintf(temp, "%s+%s", expr, temp2);
1389                                else
1390                                    sprintf(temp, "%s%s", expr, temp2);
1391                            }
1392                            assert (strlen(temp) < 1000);
1393                            setElement(iRow, j, temp);
1394                            if (amplInfo->logLevel > 1)
1395                                printf("el for row %d column c%7.7d is %s\n", iRow, j, temp);
1396                        }
1397                    }
1398                } else {
1399                    // objective
1400                    for (int j = 0; j < n_var; j++) {
1401                        for (int k = start[j]; k < start[j+1]; k++) {
1402                            int kColumn = row[k];
1403                            double value = element[k];
1404                            // ampl gives twice with assumed 0.5
1405                            if (kColumn < j)
1406                                continue;
1407                            else if (kColumn == j)
1408                                value *= 0.5;
1409                            const char * expr = getColumnObjectiveAsString(j);
1410                            double constant = 0.0;
1411                            bool linear;
1412                            if (expr && strcmp(expr, "Numeric")) {
1413                                linear = false;
1414                            } else {
1415                                constant = getColumnObjective(j);
1416                                linear = true;
1417                            }
1418                            char temp[1000];
1419                            char temp2[30];
1420                            if (value == 1.0)
1421                                sprintf(temp2, "c%7.7d", kColumn);
1422                            else
1423                                sprintf(temp2, "%g*c%7.7d", value, kColumn);
1424                            if (linear) {
1425                                if (!constant)
1426                                    strcpy(temp, temp2);
1427                                else if (value > 0.0)
1428                                    sprintf(temp, "%g+%s", constant, temp2);
1429                                else
1430                                    sprintf(temp, "%g%s", constant, temp2);
1431                            } else {
1432                                if (value > 0.0)
1433                                    sprintf(temp, "%s+%s", expr, temp2);
1434                                else
1435                                    sprintf(temp, "%s%s", expr, temp2);
1436                            }
1437                            assert (strlen(temp) < 1000);
1438                            setObjective(j, temp);
1439                            if (amplInfo->logLevel > 1)
1440                                printf("el for objective column c%7.7d is %s\n", j, temp);
1441                        }
1442                    }
1443                }
1444            }
1445        }
1446        if (nOdd) {
1447            printf("%d non-linear constraints could not be converted to quadratic\n", nOdd);
1448            exit(77);
1449        }
1450    }
1451    free(colqp);
1452    free(z);
1453    // see if any sos
1454    {
1455        char *sostype;
1456        int nsosnz, *sosbeg, *sosind, * sospri;
1457        double *sosref;
1458        int nsos;
1459        int i = ASL_suf_sos_explict_free;
1460        int copri[2], **p_sospri;
1461        copri[0] = 0;
1462        copri[1] = 0;
1463        p_sospri = &sospri;
1464        nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1465                       &sosbeg, &sosind, &sosref);
1466        if (nsos) {
1467            numberSOS_ = nsos;
1468            typeSOS_ = new int [numberSOS_];
1469            prioritySOS_ = new int [numberSOS_];
1470            startSOS_ = new int [numberSOS_+1];
1471            memberSOS_ = new int[nsosnz];
1472            referenceSOS_ = new double [nsosnz];
1473            sos_kludge(nsos, sosbeg, sosref, sosind);
1474            for (int i = 0; i < nsos; i++) {
1475                int ichar = sostype[i];
1476                assert (ichar == '1' || ichar == '2');
1477                typeSOS_[i] = ichar - '0';
1478            }
1479            memcpy(prioritySOS_, sospri, nsos*sizeof(int));
1480            memcpy(startSOS_, sosbeg, (nsos + 1)*sizeof(int));
1481            memcpy(memberSOS_, sosind, nsosnz*sizeof(int));
1482            memcpy(referenceSOS_, sosref, nsosnz*sizeof(double));
1483        }
1484    }
1485}
1486#else
1487#include "Cbc_ampl.h"
1488int
1489readAmpl(ampl_info * , int , char **, void ** )
1490{
1491    return 0;
1492}
1493void freeArrays1(ampl_info *)
1494{
1495}
1496void freeArrays2(ampl_info *)
1497{
1498}
1499void freeArgs(ampl_info * )
1500{
1501}
1502int ampl_obj_prec()
1503{
1504    return 0;
1505}
1506void writeAmpl(ampl_info * )
1507{
1508}
1509#endif
1510
Note: See TracBrowser for help on using the repository browser.