source: stable/2.8/Cbc/src/Cbc_ampl.cpp @ 2034

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

allow longer row in ampl

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 53.1 KB
Line 
1/* $Id: Cbc_ampl.cpp 2034 2014-05-13 07:25:03Z 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    int precision = obj_prec();
701    if (precision<=0)
702        precision=15;
703    return precision;
704}
705void writeAmpl(ampl_info * info)
706{
707    char buf[1000];
708    typedef struct {
709        const char *msg;
710        int code;
711        int wantObj;
712    } Sol_info;
713    static Sol_info solinfo[] = {
714        { "optimal solution",                   000, 1 },
715        { "infeasible",                         200, 1 },
716        { "unbounded",                          300, 0 },
717        { "iteration limit etc",                        400, 1 },
718        { "solution limit",                             401, 1 },
719        { "ran out of space",                   500, 0 },
720        { "status unknown",                             501, 1 },
721        { "bug!",                                       502, 0 },
722        { "best MIP solution so far restored",  101, 1 },
723        { "failed to restore best MIP solution",        503, 1 },
724        { "optimal (?) solution",                       100, 1 }
725    };
726    /* convert status - need info on map */
727    static int map[] = {0, 3, 4, 1};
728    sprintf(buf, "%s %s", Oinfo.bsname, info->buffer);
729    solve_result_num = solinfo[info->problemStatus].code;
730    if (info->columnStatus) {
731        stat_map(info->columnStatus, n_var, map, 4, "outgoing columnStatus");
732        stat_map(info->rowStatus, n_con, map, 4, "outgoing rowStatus");
733        suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
734        suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
735    }
736    write_sol(buf, info->primalSolution, info->dualSolution, &Oinfo);
737}
738/* Read a problem from AMPL nl file
739 */
740CoinModel::CoinModel( int nonLinear, const char * fileName, const void * info)
741        :  CoinBaseModel(),
742        maximumRows_(0),
743        maximumColumns_(0),
744        numberElements_(0),
745        maximumElements_(0),
746        numberQuadraticElements_(0),
747        maximumQuadraticElements_(0),
748        rowLower_(NULL),
749        rowUpper_(NULL),
750        rowType_(NULL),
751        objective_(NULL),
752        columnLower_(NULL),
753        columnUpper_(NULL),
754        integerType_(NULL),
755        columnType_(NULL),
756        start_(NULL),
757        elements_(NULL),
758        packedMatrix_(NULL),
759        quadraticElements_(NULL),
760        sortIndices_(NULL),
761        sortElements_(NULL),
762        sortSize_(0),
763        sizeAssociated_(0),
764        associated_(NULL),
765        numberSOS_(0),
766        startSOS_(NULL),
767        memberSOS_(NULL),
768        typeSOS_(NULL),
769        prioritySOS_(NULL),
770        referenceSOS_(NULL),
771        priority_(NULL),
772        cut_(NULL),
773        moreInfo_(NULL),
774        type_(-1),
775        noNames_(false),
776        links_(0)
777{
778    problemName_ = "";
779    int status = 0;
780    if (!strcmp(fileName, "-") || !strcmp(fileName, "stdin")) {
781        // stdin
782    } else {
783        std::string name = fileName;
784        bool readable = fileCoinReadable(name);
785        if (!readable) {
786            std::cerr << "Unable to open file "
787                      << fileName << std::endl;
788            status = -1;
789        }
790    }
791    if (!status) {
792        gdb(nonLinear, fileName, info);
793    }
794}
795#ifdef JJF_ZERO
796static real
797qterm(ASL *asl, fint *colq, fint *rowq, real *delsq)
798{
799    double t, t1, *x, *x0, *xe;
800    fint *rq0, *rqe;
801
802    t = 0.;
803    x0 = x = X0;
804    xe = x + n_var;
805    rq0 = rowq;
806    while (x < xe) {
807        t1 = *x++;
808        rqe = rq0 + *++colq;
809        while (rowq < rqe)
810            t += t1 * x0[*rowq++]**delsq++;
811    }
812    return 0.5 * t;
813}
814#endif
815// stolen from IPopt with changes
816typedef struct {
817    double obj_sign_;
818    ASL_pfgh * asl_;
819    double * non_const_x_;
820    int * column_; // for jacobian
821    int * rowStart_;
822    double * gradient_;
823    double * constraintValues_;
824    int nz_h_full_; // number of nonzeros in hessian
825    int nerror_;
826    bool objval_called_with_current_x_;
827    bool conval_called_with_current_x_;
828    bool jacval_called_with_current_x_;
829} CbcAmplInfo;
830
831void
832CoinModel::gdb( int nonLinear, const char * fileName, const void * info)
833{
834    const ampl_info * amplInfo = (const ampl_info *) info;
835    ograd *og = NULL;
836    int i;
837    SufDesc *csd = NULL;
838    SufDesc *rsd = NULL;
839    /*bool *basis, *lower;*/
840    /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
841    //char tempBuffer[20];
842    double * objective = NULL;
843    double * columnLower = NULL;
844    double * columnUpper = NULL;
845    double * rowLower = NULL;
846    double * rowUpper = NULL;
847    int * columnStatus = NULL;
848    int * rowStatus = NULL;
849    int numberRows = -1;
850    int numberColumns = -1;
851    int numberElements = -1;
852    int numberBinary = -1;
853    int numberIntegers = -1;
854    int numberAllNonLinearBoth = 0;
855    int numberIntegerNonLinearBoth = 0;
856    int numberAllNonLinearConstraints = 0;
857    int numberIntegerNonLinearConstraints = 0;
858    int numberAllNonLinearObjective = 0;
859    int numberIntegerNonLinearObjective = 0;
860    double * primalSolution = NULL;
861    double direction = 1.0;
862    char * stub = strdup(fileName);
863    CoinPackedMatrix matrixByRow;
864    fint ** colqp = NULL;
865    int *z = NULL;
866    if (nonLinear == 0) {
867        // linear
868        asl = ASL_alloc(ASL_read_f);
869        nl = jac0dim(stub, 0);
870        free(stub);
871        suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
872
873        /* set A_vals to get the constraints column-wise (malloc so can be freed) */
874        A_vals = (double *) malloc(nzc * sizeof(double));
875        if (!A_vals) {
876            printf("no memory\n");
877            return ;
878        }
879        /* say we want primal solution */
880        want_xpi0 = 1;
881        /* for basis info */
882        columnStatus = (int *) malloc(n_var * sizeof(int));
883        rowStatus = (int *) malloc(n_con * sizeof(int));
884        csd = suf_iput("sstatus", ASL_Sufkind_var, columnStatus);
885        rsd = suf_iput("sstatus", ASL_Sufkind_con, rowStatus);
886        /* read linear model*/
887        f_read(nl, 0);
888        // see if any sos
889        if (true) {
890            char *sostype;
891            int nsosnz, *sosbeg, *sosind, * sospri;
892            double *sosref;
893            int nsos;
894            int i = ASL_suf_sos_explict_free;
895            int copri[2], **p_sospri;
896            copri[0] = 0;
897            copri[1] = 0;
898            p_sospri = &sospri;
899            nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
900                           &sosbeg, &sosind, &sosref);
901            if (nsos) {
902                abort();
903#ifdef JJF_ZERO
904                info->numberSos = nsos;
905                info->sosType = (char *) malloc(nsos);
906                info->sosPriority = (int *) malloc(nsos * sizeof(int));
907                info->sosStart = (int *) malloc((nsos + 1) * sizeof(int));
908                info->sosIndices = (int *) malloc(nsosnz * sizeof(int));
909                info->sosReference = (double *) malloc(nsosnz * sizeof(double));
910                sos_kludge(nsos, sosbeg, sosref, sosind);
911                for (int i = 0; i < nsos; i++) {
912                    int ichar = sostype[i];
913                    assert (ichar == '1' || ichar == '2');
914                    info->sosType[i] = ichar - '0';
915                }
916                memcpy(info->sosPriority, sospri, nsos*sizeof(int));
917                memcpy(info->sosStart, sosbeg, (nsos + 1)*sizeof(int));
918                memcpy(info->sosIndices, sosind, nsosnz*sizeof(int));
919                memcpy(info->sosReference, sosref, nsosnz*sizeof(double));
920#endif
921            }
922        }
923
924        /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
925        //Oinfo.uinfo = tempBuffer;
926        //if (getopts(argv, &Oinfo))
927        //return 1;
928        /* objective*/
929        objective = (double *) malloc(n_var * sizeof(double));
930        for (i = 0; i < n_var; i++)
931            objective[i] = 0.0;
932        if (n_obj) {
933            for (og = Ograd[0]; og; og = og->next)
934                objective[og->varno] = og->coef;
935        }
936        if (objtype[0])
937            direction = -1.0;
938        else
939            direction = 1.0;
940        objectiveOffset_ = objconst(0);
941        /* Column bounds*/
942        columnLower = (double *) malloc(n_var * sizeof(double));
943        columnUpper = (double *) malloc(n_var * sizeof(double));
944        for (i = 0; i < n_var; i++) {
945            columnLower[i] = LUv[2*i];
946            if (columnLower[i] <= negInfinity)
947                columnLower[i] = -COIN_DBL_MAX;
948            columnUpper[i] = LUv[2*i+1];
949            if (columnUpper[i] >= Infinity)
950                columnUpper[i] = COIN_DBL_MAX;
951        }
952        /* Row bounds*/
953        rowLower = (double *) malloc(n_con * sizeof(double));
954        rowUpper = (double *) malloc(n_con * sizeof(double));
955        for (i = 0; i < n_con; i++) {
956            rowLower[i] = LUrhs[2*i];
957            if (rowLower[i] <= negInfinity)
958                rowLower[i] = -COIN_DBL_MAX;
959            rowUpper[i] = LUrhs[2*i+1];
960            if (rowUpper[i] >= Infinity)
961                rowUpper[i] = COIN_DBL_MAX;
962        }
963        numberRows = n_con;
964        numberColumns = n_var;
965        numberElements = nzc;
966        numberBinary = nbv;
967        numberIntegers = niv;
968        /* put in primalSolution if exists */
969        if (X0) {
970            primalSolution = (double *) malloc(n_var * sizeof(double));
971            memcpy( primalSolution, X0, n_var*sizeof(double));
972        }
973        //double * dualSolution=NULL;
974        if (niv + nbv > 0)
975            mip_stuff(); // get any extra info
976        if ((!(niv + nbv) && (csd->kind & ASL_Sufkind_input))
977                || (rsd->kind & ASL_Sufkind_input)) {
978            /* convert status - need info on map */
979            static int map[] = {1, 3, 1, 1, 2, 1, 1};
980            stat_map(columnStatus, n_var, map, 6, "incoming columnStatus");
981            stat_map(rowStatus, n_con, map, 6, "incoming rowStatus");
982        } else {
983            /* all slack basis */
984            // leave status for output */
985#ifdef JJF_ZERO
986            free(rowStatus);
987            rowStatus = NULL;
988            free(columnStatus);
989            columnStatus = NULL;
990#endif
991        }
992        CoinPackedMatrix columnCopy(true, numberRows, numberColumns, numberElements,
993                                    A_vals, A_rownos, A_colstarts, NULL);
994        matrixByRow.reverseOrderedCopyOf(columnCopy);
995    } else if (nonLinear == 1) {
996        // quadratic
997        asl = ASL_alloc(ASL_read_fg);
998        nl = jac0dim(stub, (ftnlen) strlen(stub));
999        free(stub);
1000        suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
1001        /* read  model*/
1002        X0 = (double*) malloc(n_var * sizeof(double));
1003        CoinZeroN(X0, n_var);
1004        qp_read(nl, 0);
1005        assert (n_obj == 1);
1006        int nz = 1 + n_con;
1007        colqp = (fint**) malloc(nz * (2 * sizeof(int*)
1008                                      + sizeof(double*)));
1009        fint ** rowqp = colqp + nz;
1010        double ** delsqp = (double **)(rowqp + nz);
1011        z = (int*) malloc(nz * sizeof(int));
1012        for (i = 0; i <= n_con; i++) {
1013            z[i] = nqpcheck(-i, rowqp + i, colqp + i, delsqp + i);
1014        }
1015        qp_opify();
1016        /* objective*/
1017        objective = (double *) malloc(n_var * sizeof(double));
1018        for (i = 0; i < n_var; i++)
1019            objective[i] = 0.0;
1020        if (n_obj) {
1021            for (og = Ograd[0]; og; og = og->next)
1022                objective[og->varno] = og->coef;
1023        }
1024        if (objtype[0])
1025            direction = -1.0;
1026        else
1027            direction = 1.0;
1028        objectiveOffset_ = objconst(0);
1029        /* Column bounds*/
1030        columnLower = (double *) malloc(n_var * sizeof(double));
1031        columnUpper = (double *) malloc(n_var * sizeof(double));
1032        for (i = 0; i < n_var; i++) {
1033            columnLower[i] = LUv[2*i];
1034            if (columnLower[i] <= negInfinity)
1035                columnLower[i] = -COIN_DBL_MAX;
1036            columnUpper[i] = LUv[2*i+1];
1037            if (columnUpper[i] >= Infinity)
1038                columnUpper[i] = COIN_DBL_MAX;
1039        }
1040        // Build by row from scratch
1041        //matrixByRow.reserve(n_var,nzc,true);
1042        // say row orderded
1043        matrixByRow.transpose();
1044        /* Row bounds*/
1045        rowLower = (double *) malloc(n_con * sizeof(double));
1046        rowUpper = (double *) malloc(n_con * sizeof(double));
1047        CoinBigIndex * rowStart = new CoinBigIndex [n_con+1];
1048        int * column = new int [nzc];
1049        double * element = new double [nzc];
1050        rowStart[0] = 0;
1051        numberElements = 0;
1052        for (i = 0; i < n_con; i++) {
1053            rowLower[i] = LUrhs[2*i];
1054            if (rowLower[i] <= negInfinity)
1055                rowLower[i] = -COIN_DBL_MAX;
1056            rowUpper[i] = LUrhs[2*i+1];
1057            if (rowUpper[i] >= Infinity)
1058                rowUpper[i] = COIN_DBL_MAX;
1059            for (cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
1060                column[numberElements] = cg->varno;
1061                element[numberElements++] = cg->coef;
1062            }
1063            rowStart[i+1] = numberElements;
1064        }
1065        assert (numberElements == nzc);
1066        matrixByRow.appendRows(n_con, rowStart, column, element);
1067        delete [] rowStart;
1068        delete [] column;
1069        delete [] element;
1070        numberRows = n_con;
1071        numberColumns = n_var;
1072        //numberElements=nzc;
1073        numberBinary = nbv;
1074        numberIntegers = niv;
1075        numberAllNonLinearBoth = nlvb;
1076        numberIntegerNonLinearBoth = nlvbi;
1077        numberAllNonLinearConstraints = nlvc;
1078        numberIntegerNonLinearConstraints = nlvci;
1079        numberAllNonLinearObjective = nlvo;
1080        numberIntegerNonLinearObjective = nlvoi;
1081        /* say we want primal solution */
1082        want_xpi0 = 1;
1083        //double * dualSolution=NULL;
1084        // save asl
1085        // Fix memory leak one day
1086        CbcAmplInfo * info = new CbcAmplInfo;
1087        //amplGamsData_ = info;
1088        info->asl_ = NULL; // as wrong form asl;
1089        info->nz_h_full_ = -1; // number of nonzeros in hessian
1090        info->objval_called_with_current_x_ = false;
1091        info->nerror_ = 0;
1092        info->obj_sign_ = direction;
1093        info->conval_called_with_current_x_ = false;
1094        info->non_const_x_ = NULL;
1095        info->jacval_called_with_current_x_ = false;
1096        info->rowStart_ = NULL;
1097        info->column_ = NULL;
1098        info->gradient_ = NULL;
1099        info->constraintValues_ = NULL;
1100    } else if (nonLinear == 2) {
1101        // General nonlinear!
1102        //ASL_pfgh* asl = (ASL_pfgh*)ASL_alloc(ASL_read_pfgh);
1103        asl = ASL_alloc(ASL_read_pfgh);
1104        nl = jac0dim(stub, (ftnlen) strlen(stub));
1105        free(stub);
1106        suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
1107        /* read  model*/
1108        X0 = (double*) malloc(n_var * sizeof(double));
1109        CoinZeroN(X0, n_var);
1110        // code stolen from Ipopt
1111        int retcode = pfgh_read(nl, ASL_return_read_err | ASL_findgroups);
1112
1113        switch (retcode) {
1114        case ASL_readerr_none : {}
1115        break;
1116        case ASL_readerr_nofile : {
1117            printf( "Cannot open .nl file\n");
1118            exit(-1);
1119        }
1120        break;
1121        case ASL_readerr_nonlin : {
1122            assert(false); // this better not be an error!
1123            printf( "model involves nonlinearities (ed0read)\n");
1124            exit(-1);
1125        }
1126        break;
1127        case  ASL_readerr_argerr : {
1128            printf( "user-defined function with bad args\n");
1129            exit(-1);
1130        }
1131        break;
1132        case ASL_readerr_unavail : {
1133            printf( "user-defined function not available\n");
1134            exit(-1);
1135        }
1136        break;
1137        case ASL_readerr_corrupt : {
1138            printf( "corrupt .nl file\n");
1139            exit(-1);
1140        }
1141        break;
1142        case ASL_readerr_bug : {
1143            printf( "bug in .nl reader\n");
1144            exit(-1);
1145        }
1146        break;
1147        case ASL_readerr_CLP : {
1148            printf( "ASL error message: \"solver cannot handle CLP extensions\"\n");
1149            exit(-1);
1150        }
1151        break;
1152        default: {
1153            printf( "Unknown error in stub file read. retcode = %d\n", retcode);
1154            exit(-1);
1155        }
1156        break;
1157        }
1158
1159        // see "changes" in solvers directory of ampl code...
1160        hesset(1, 0, 1, 0, nlc);
1161
1162        assert (n_obj == 1);
1163        // find the nonzero structure for the hessian
1164        // parameters to sphsetup:
1165        int coeff_obj = 1; // coefficient of the objective fn ???
1166        int mult_supplied = 1; // multipliers will be supplied
1167        int uptri = 1; // only need the upper triangular part
1168        // save asl
1169        // Fix memory leak one day
1170        CbcAmplInfo * info = new CbcAmplInfo;
1171        moreInfo_ = (void *) info;
1172        //amplGamsData_ = info;
1173        info->asl_ = (ASL_pfgh *) asl;
1174        // This is not easy to get from ampl so save
1175        info->nz_h_full_ = sphsetup(-1, coeff_obj, mult_supplied, uptri);
1176        info->objval_called_with_current_x_ = false;
1177        info->nerror_ = 0;
1178        info->obj_sign_ = direction;
1179        info->conval_called_with_current_x_ = false;
1180        info->non_const_x_ = NULL;
1181        info->jacval_called_with_current_x_ = false;
1182        // Look at nonlinear
1183        if (nzc) {
1184            n_conjac[1] = nlc; // just nonlinear
1185            int * rowStart = new int [nlc+1];
1186            info->rowStart_ = rowStart;
1187            // See how many
1188            int  current_nz = 0;
1189            for (int i = 0; i < nlc; i++) {
1190                for (cgrad* cg = Cgrad[i]; cg; cg = cg->next) {
1191                    current_nz++;
1192                }
1193            }
1194            // setup the structure
1195            int * column = new int [current_nz];
1196            info->column_ = column;
1197            current_nz = 0;
1198            rowStart[0] = 0;
1199            for (int i = 0; i < nlc; i++) {
1200                for (cgrad* cg = Cgrad[i]; cg; cg = cg->next) {
1201                    cg->goff = current_nz;
1202                    //iRow[cg->goff] = i ;
1203                    //jCol[cg->goff] = cg->varno + 1;
1204                    column[cg->goff] = cg->varno ;
1205                    current_nz++;
1206                }
1207                rowStart[i+1] = current_nz;
1208            }
1209            info->gradient_ = new double [nzc];
1210            info->constraintValues_ = new double [nlc];
1211        }
1212        /* objective*/
1213        objective = (double *) malloc(n_var * sizeof(double));
1214        for (i = 0; i < n_var; i++)
1215            objective[i] = 0.0;
1216        if (n_obj) {
1217            for (og = Ograd[0]; og; og = og->next)
1218                objective[og->varno] = og->coef;
1219        }
1220        if (objtype[0])
1221            direction = -1.0;
1222        else
1223            direction = 1.0;
1224        objectiveOffset_ = objconst(0);
1225        /* Column bounds*/
1226        columnLower = (double *) malloc(n_var * sizeof(double));
1227        columnUpper = (double *) malloc(n_var * sizeof(double));
1228        for (i = 0; i < n_var; i++) {
1229            columnLower[i] = LUv[2*i];
1230            if (columnLower[i] <= negInfinity)
1231                columnLower[i] = -COIN_DBL_MAX;
1232            columnUpper[i] = LUv[2*i+1];
1233            if (columnUpper[i] >= Infinity)
1234                columnUpper[i] = COIN_DBL_MAX;
1235        }
1236        // Build by row from scratch
1237        //matrixByRow.reserve(n_var,nzc,true);
1238        // say row orderded
1239        matrixByRow.transpose();
1240        CoinBigIndex * rowStart = new CoinBigIndex [n_con+1];
1241        int * column = new int [nzc];
1242        double * element = new double [nzc];
1243        rowStart[0] = 0;
1244        numberElements = 0;
1245        /* Row bounds*/
1246        rowLower = (double *) malloc(n_con * sizeof(double));
1247        rowUpper = (double *) malloc(n_con * sizeof(double));
1248        for (i = 0; i < n_con; i++) {
1249            rowLower[i] = LUrhs[2*i];
1250            if (rowLower[i] <= negInfinity)
1251                rowLower[i] = -COIN_DBL_MAX;
1252            rowUpper[i] = LUrhs[2*i+1];
1253            if (rowUpper[i] >= Infinity)
1254                rowUpper[i] = COIN_DBL_MAX;
1255            for (cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
1256                column[numberElements] = cg->varno;
1257                double value = cg->coef;
1258                if (!value)
1259                    value = -1.2345e-29;
1260                element[numberElements++] = value;
1261            }
1262            rowStart[i+1] = numberElements;
1263        }
1264        assert (numberElements == nzc);
1265        matrixByRow.appendRows(n_con, rowStart, column, element);
1266        delete [] rowStart;
1267        delete [] column;
1268        delete [] element;
1269        numberRows = n_con;
1270        numberColumns = n_var;
1271        numberElements = nzc;
1272        numberBinary = nbv;
1273        numberIntegers = niv;
1274        numberAllNonLinearBoth = nlvb;
1275        numberIntegerNonLinearBoth = nlvbi;
1276        numberAllNonLinearConstraints = nlvc;
1277        numberIntegerNonLinearConstraints = nlvci;
1278        numberAllNonLinearObjective = nlvo;
1279        numberIntegerNonLinearObjective = nlvoi;
1280        /* say we want primal solution */
1281        want_xpi0 = 1;
1282        //double * dualSolution=NULL;
1283    } else {
1284        abort();
1285    }
1286    // set problem name
1287    problemName_ = "???";
1288
1289    // Build by row from scratch
1290    const double * element = matrixByRow.getElements();
1291    const int * column = matrixByRow.getIndices();
1292    const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
1293    const int * rowLength = matrixByRow.getVectorLengths();
1294    for (i = 0; i < numberRows; i++) {
1295        addRow(rowLength[i], column + rowStart[i],
1296               element + rowStart[i], rowLower[i], rowUpper[i]);
1297    }
1298    // Now do column part
1299    for (i = 0; i < numberColumns; i++) {
1300        setColumnBounds(i, columnLower[i], columnUpper[i]);
1301        setColumnObjective(i, objective[i]);
1302    }
1303    for ( i = numberColumns - numberBinary - numberIntegers;
1304            i < numberColumns; i++) {
1305        setColumnIsInteger(i, true);
1306    }
1307    // and non linear
1308    for (i = numberAllNonLinearBoth - numberIntegerNonLinearBoth;
1309            i < numberAllNonLinearBoth; i++) {
1310        setColumnIsInteger(i, true);
1311    }
1312    for (i = numberAllNonLinearConstraints - numberIntegerNonLinearConstraints;
1313            i < numberAllNonLinearConstraints; i++) {
1314        setColumnIsInteger(i, true);
1315    }
1316    for (i = numberAllNonLinearObjective - numberIntegerNonLinearObjective;
1317            i < numberAllNonLinearObjective; i++) {
1318        setColumnIsInteger(i, true);
1319    }
1320    free(columnLower);
1321    free(columnUpper);
1322    free(rowLower);
1323    free(rowUpper);
1324    free(objective);
1325    // space for building a row
1326    char * temp = new char [30*numberColumns_];
1327    // do names
1328    int iRow;
1329    for (iRow = 0; iRow < numberRows_; iRow++) {
1330        char name[9];
1331        sprintf(name, "r%7.7d", iRow);
1332        setRowName(iRow, name);
1333    }
1334    int iColumn;
1335    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1336        char name[9];
1337        sprintf(name, "c%7.7d", iColumn);
1338        setColumnName(iColumn, name);
1339    }
1340    if (colqp) {
1341        // add in quadratic
1342        int nz = 1 + n_con;
1343        int nOdd = 0;
1344        fint ** rowqp = colqp + nz;
1345        double ** delsqp = (double **)(rowqp + nz);
1346        for (i = 0; i <= n_con; i++) {
1347            int nels = z[i];
1348            if (nels) {
1349                double * element = delsqp[i];
1350                int * start = (int *) colqp[i];
1351                int * row = (int *) rowqp[i];
1352                if (!element) {
1353                    // odd row - probably not quadratic
1354                    nOdd++;
1355                    continue;
1356                }
1357#ifdef JJF_ZERO
1358                printf("%d quadratic els\n", nels);
1359                for (int j = 0; j < n_var; j++) {
1360                    for (int k = start[j]; k < start[j+1]; k++)
1361                        printf("%d %d %g\n", j, row[k], element[k]);
1362                }
1363#endif
1364                if (i) {
1365                    int iRow = i - 1;
1366                    for (int j = 0; j < n_var; j++) {
1367                        for (int k = start[j]; k < start[j+1]; k++) {
1368                            int kColumn = row[k];
1369                            double value = element[k];
1370                            // ampl gives twice with assumed 0.5
1371                            if (kColumn < j)
1372                                continue;
1373                            else if (kColumn == j)
1374                                value *= 0.5;
1375                            const char * expr = getElementAsString(iRow, j);
1376                            double constant = 0.0;
1377                            bool linear;
1378                            if (expr && strcmp(expr, "Numeric")) {
1379                                linear = false;
1380                            } else {
1381                                constant = getElement(iRow, j);
1382                                linear = true;
1383                            }
1384                            char temp2[30];
1385                            if (value == 1.0)
1386                                sprintf(temp2, "c%7.7d", kColumn);
1387                            else
1388                                sprintf(temp2, "%g*c%7.7d", value, kColumn);
1389                            if (linear) {
1390                                if (!constant)
1391                                    strcpy(temp, temp2);
1392                                else if (value > 0.0)
1393                                    sprintf(temp, "%g+%s", constant, temp2);
1394                                else
1395                                    sprintf(temp, "%g%s", constant, temp2);
1396                            } else {
1397                                if (value > 0.0)
1398                                    sprintf(temp, "%s+%s", expr, temp2);
1399                                else
1400                                    sprintf(temp, "%s%s", expr, temp2);
1401                            }
1402                            assert (static_cast<int>(strlen(temp)) < 30*numberColumns_);
1403                            setElement(iRow, j, temp);
1404                            if (amplInfo->logLevel > 1)
1405                                printf("el for row %d column c%7.7d is %s\n", iRow, j, temp);
1406                        }
1407                    }
1408                } else {
1409                    // objective
1410                    for (int j = 0; j < n_var; j++) {
1411                        for (int k = start[j]; k < start[j+1]; k++) {
1412                            int kColumn = row[k];
1413                            double value = element[k];
1414                            // ampl gives twice with assumed 0.5
1415                            if (kColumn < j)
1416                                continue;
1417                            else if (kColumn == j)
1418                                value *= 0.5;
1419                            const char * expr = getColumnObjectiveAsString(j);
1420                            double constant = 0.0;
1421                            bool linear;
1422                            if (expr && strcmp(expr, "Numeric")) {
1423                                linear = false;
1424                            } else {
1425                                constant = getColumnObjective(j);
1426                                linear = true;
1427                            }
1428                            char temp2[30];
1429                            if (value == 1.0)
1430                                sprintf(temp2, "c%7.7d", kColumn);
1431                            else
1432                                sprintf(temp2, "%g*c%7.7d", value, kColumn);
1433                            if (linear) {
1434                                if (!constant)
1435                                    strcpy(temp, temp2);
1436                                else if (value > 0.0)
1437                                    sprintf(temp, "%g+%s", constant, temp2);
1438                                else
1439                                    sprintf(temp, "%g%s", constant, temp2);
1440                            } else {
1441                                if (value > 0.0)
1442                                    sprintf(temp, "%s+%s", expr, temp2);
1443                                else
1444                                    sprintf(temp, "%s%s", expr, temp2);
1445                            }
1446                            assert (static_cast<int>(strlen(temp)) < 30*numberColumns_);
1447                            setObjective(j, temp);
1448                            if (amplInfo->logLevel > 1)
1449                                printf("el for objective column c%7.7d is %s\n", j, temp);
1450                        }
1451                    }
1452                }
1453            }
1454        }
1455        if (nOdd) {
1456            printf("%d non-linear constraints could not be converted to quadratic\n", nOdd);
1457            exit(77);
1458        }
1459    }
1460    delete [] temp;
1461    free(colqp);
1462    free(z);
1463    // see if any sos
1464    {
1465        char *sostype;
1466        int nsosnz, *sosbeg, *sosind, * sospri;
1467        double *sosref;
1468        int nsos;
1469        int i = ASL_suf_sos_explict_free;
1470        int copri[2], **p_sospri;
1471        copri[0] = 0;
1472        copri[1] = 0;
1473        p_sospri = &sospri;
1474        nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1475                       &sosbeg, &sosind, &sosref);
1476        if (nsos) {
1477            numberSOS_ = nsos;
1478            typeSOS_ = new int [numberSOS_];
1479            prioritySOS_ = new int [numberSOS_];
1480            startSOS_ = new int [numberSOS_+1];
1481            memberSOS_ = new int[nsosnz];
1482            referenceSOS_ = new double [nsosnz];
1483            sos_kludge(nsos, sosbeg, sosref, sosind);
1484            for (int i = 0; i < nsos; i++) {
1485                int ichar = sostype[i];
1486                assert (ichar == '1' || ichar == '2');
1487                typeSOS_[i] = ichar - '0';
1488            }
1489            memcpy(prioritySOS_, sospri, nsos*sizeof(int));
1490            memcpy(startSOS_, sosbeg, (nsos + 1)*sizeof(int));
1491            memcpy(memberSOS_, sosind, nsosnz*sizeof(int));
1492            memcpy(referenceSOS_, sosref, nsosnz*sizeof(double));
1493        }
1494    }
1495}
1496#else
1497#include "Cbc_ampl.h"
1498int
1499readAmpl(ampl_info * , int , char **, void ** )
1500{
1501    return 0;
1502}
1503void freeArrays1(ampl_info *)
1504{
1505}
1506void freeArrays2(ampl_info *)
1507{
1508}
1509void freeArgs(ampl_info * )
1510{
1511}
1512int ampl_obj_prec()
1513{
1514    return 0;
1515}
1516void writeAmpl(ampl_info * )
1517{
1518}
1519#endif
1520
Note: See TracBrowser for help on using the repository browser.