source: branches/sandbox/Cbc/src/Cbc_ampl.cpp @ 1289

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

Remove extra Bjarni

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