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

Last change on this file since 1286 was 1286, checked in by EdwinStraver, 10 years ago

Changed formatting using AStyle -A4 -p

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