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

Last change on this file since 2345 was 2345, checked in by forrest, 3 years ago

ampl with COIN_BIG_INDEX

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