source: releases/2.8.9/Cbc/src/Cbc_ampl.cpp @ 2028

Last change on this file since 2028 was 2028, checked in by forrest, 5 years ago

fix ampl uninitialized memory

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