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

Last change on this file since 1873 was 1873, checked in by stefan, 7 years ago

sync with trunk rev 1872

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