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

Last change on this file since 1926 was 1926, checked in by stefan, 5 years ago

sync with trunk rev 1925

  • 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 1926 2013-05-24 10:19:49Z 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    20130502
159};
160// strdup used to avoid g++ compiler warning
161static SufDecl suftab[] = {
162#ifdef JJF_ZERO
163    { const_cast<char*>("current"), 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
164    { const_cast<char*>("current"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
165    { const_cast<char*>("direction"), 0, ASL_Sufkind_var },
166    { const_cast<char*>("down"), 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
167    { const_cast<char*>("down"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
168    { const_cast<char*>("priority"), 0, ASL_Sufkind_var },
169#endif
170    { const_cast<char*>("cut"), 0, ASL_Sufkind_con },
171    { const_cast<char*>("direction"), 0, ASL_Sufkind_var },
172    { const_cast<char*>("downPseudocost"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
173    { const_cast<char*>("priority"), 0, ASL_Sufkind_var },
174    { const_cast<char*>("ref"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
175    { const_cast<char*>("sos"), 0, ASL_Sufkind_var },
176    { const_cast<char*>("sos"), 0, ASL_Sufkind_con },
177    { const_cast<char*>("sosno"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
178    { const_cast<char*>("sosref"), 0, ASL_Sufkind_var | ASL_Sufkind_real },
179    { const_cast<char*>("special"), 0, ASL_Sufkind_var },
180    { const_cast<char*>("special"), 0, ASL_Sufkind_con },
181    /*{ const_cast<char*>("special"), 0, ASL_Sufkind_con },*/
182    { const_cast<char*>("sstatus"), 0, ASL_Sufkind_var, 0 },
183    { const_cast<char*>("sstatus"), 0, ASL_Sufkind_con, 0 },
184    { const_cast<char*>("upPseudocost"), 0, ASL_Sufkind_var | ASL_Sufkind_real }
185#ifdef JJF_ZERO
186    { const_cast<char*>("unbdd"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly},
187    { const_cast<char*>("up"), 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
188    { const_cast<char*>("up"), 0, ASL_Sufkind_var | ASL_Sufkind_outonly }
189#endif
190};
191#include "float.h"
192#include "limits.h"
193static ASL *asl = NULL;
194static FILE *nl = NULL;
195
196static void
197mip_stuff(void)
198{
199    int i;
200    double *pseudoUp, *pseudoDown;
201    int *priority, *direction;
202    // To label cuts (there will be other uses for special)
203    int *cut;
204    // To label special variables - at present 1= must be >= 1 or <= -1
205    int * special;
206    SufDesc *dpup, *dpdown, *dpri, *ddir, *dcut, *dspecial;
207
208    ddir = suf_get("direction", ASL_Sufkind_var);
209    direction = ddir->u.i;
210    dpri = suf_get("priority", ASL_Sufkind_var);
211    priority = dpri->u.i;
212    dspecial = suf_get("special", ASL_Sufkind_con);
213    dcut = suf_get("cut", ASL_Sufkind_con);
214    cut = dcut->u.i;
215    if (!cut) {
216        // try special
217        dcut = suf_get("special", ASL_Sufkind_con);
218        cut = dcut->u.i;
219    }
220    dspecial = suf_get("special", ASL_Sufkind_var);
221    special = dspecial->u.i;
222    dpdown = suf_get("downPseudocost", ASL_Sufkind_var);
223    pseudoDown = dpdown->u.r;
224    dpup = suf_get("upPseudocost", ASL_Sufkind_var);
225    pseudoUp = dpup->u.r;
226    assert(saveInfo);
227    int numberColumns = saveInfo->numberColumns;
228    if (direction) {
229        int baddir = 0;
230        saveInfo->branchDirection = (int *) malloc(numberColumns * sizeof(int));
231        for (i = 0; i < numberColumns; i++) {
232            int value = direction[i];
233            if (value < -1 || value > 1) {
234                baddir++;
235                value = 0;
236            }
237            saveInfo->branchDirection[i] = value;
238        }
239        if (baddir)
240            fprintf(Stderr,
241                    "Treating %d .direction values outside [-1, 1] as 0.\n",
242                    baddir);
243    }
244    if (priority) {
245        int badpri = 0;
246        saveInfo->priorities = (int *) malloc(numberColumns * sizeof(int));
247        for (i = 0; i < numberColumns; i++) {
248            int value = priority[i];
249            if (value < 0) {
250                badpri++;
251                value = 0;
252            }
253            saveInfo->priorities[i] = value;
254        }
255        if (badpri)
256            fprintf(Stderr,
257                    "Treating %d negative .priority values as 0\n",
258                    badpri);
259    }
260    if (special) {
261        int badspecial = 0;
262        saveInfo->special = (int *) malloc(numberColumns * sizeof(int));
263        for (i = 0; i < numberColumns; i++) {
264            int value = special[i];
265            if (value < 0) {
266                badspecial++;
267                value = 0;
268            }
269            saveInfo->special[i] = value;
270        }
271        if (badspecial)
272            fprintf(Stderr,
273                    "Treating %d negative special values as 0\n",
274                    badspecial);
275    }
276    int numberRows = saveInfo->numberRows;
277    if (cut) {
278        int badcut = 0;
279        saveInfo->cut = (int *) malloc(numberRows * sizeof(int));
280        for (i = 0; i < numberRows; i++) {
281            int value = cut[i];
282            if (value < 0) {
283                badcut++;
284                value = 0;
285            }
286            saveInfo->cut[i] = value;
287        }
288        if (badcut)
289            fprintf(Stderr,
290                    "Treating %d negative cut values as 0\n",
291                    badcut);
292    }
293    if (pseudoDown || pseudoUp) {
294        int badpseudo = 0;
295        if (!pseudoDown || !pseudoUp)
296            fprintf(Stderr,
297                    "Only one set of pseudocosts - assumed same\n");
298        saveInfo->pseudoDown = (double *) malloc(numberColumns * sizeof(double));
299        saveInfo->pseudoUp = (double *) malloc(numberColumns * sizeof(double));
300        for (i = 0; i < numberColumns; i++) {
301            double valueD = 0.0, valueU = 0.0;
302            if (pseudoDown) {
303                valueD = pseudoDown[i];
304                if (valueD < 0) {
305                    badpseudo++;
306                    valueD = 0.0;
307                }
308            }
309            if (pseudoUp) {
310                valueU = pseudoUp[i];
311                if (valueU < 0) {
312                    badpseudo++;
313                    valueU = 0.0;
314                }
315            }
316            if (!valueD)
317                valueD = valueU;
318            if (!valueU)
319                valueU = valueD;
320            saveInfo->pseudoDown[i] = valueD;
321            saveInfo->pseudoUp[i] = valueU;
322        }
323        if (badpseudo)
324            fprintf(Stderr,
325                    "Treating %d negative pseudoCosts as 0.0\n", badpseudo);
326    }
327}
328static void
329stat_map(int *stat, int n, int *map, int mx, const char *what)
330{
331    int bad, i, i1 = 0, j, j1 = 0;
332    static char badfmt[] = "Coin driver: %s[%d] = %d\n";
333
334    for (i = bad = 0; i < n; i++) {
335        if ((j = stat[i]) >= 0 && j <= mx)
336            stat[i] = map[j];
337        else {
338            stat[i] = 0;
339            i1 = i;
340            j1 = j;
341            if (!bad++)
342                fprintf(Stderr, badfmt, what, i, j);
343        }
344    }
345    if (bad > 1) {
346        if (bad == 2)
347            fprintf(Stderr, badfmt, what, i1, j1);
348        else
349            fprintf(Stderr,
350                    "Coin driver: %d messages about bad %s values suppressed.\n",
351                    bad - 1, what);
352    }
353}
354
355int
356readAmpl(ampl_info * info, int argc, char **argv, void ** coinModel)
357{
358    char *stub;
359    ograd *og;
360    int i;
361    SufDesc *csd;
362    SufDesc *rsd;
363    /*bool *basis, *lower;*/
364    /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
365    char * environment = getenv("cbc_options");
366    char tempBuffer[20];
367    double * obj;
368    double * columnLower;
369    double * columnUpper;
370    double * rowLower;
371    double * rowUpper;
372    char ** saveArgv = argv;
373    char fileName[1000];
374    if (argc > 1)
375        strcpy(fileName, argv[1]);
376    else
377        fileName[0] = '\0';
378    int nonLinearType = -1;
379    // testosi parameter - if >= 10 then go in through coinModel
380    for (i = 1; i < argc; i++) {
381        if (!strncmp(argv[i], "testosi", 7)) {
382            char * equals = strchr(argv[i], '=');
383            if (equals && atoi(equals + 1) >= 10 && atoi(equals + 1) <= 20) {
384                nonLinearType = atoi(equals + 1);
385                break;
386            }
387        }
388    }
389    int saveArgc = argc;
390    if (info->numberRows != -1234567)
391        memset(info, 0, sizeof(ampl_info)); // overwrite unless magic number set
392    /* save so can be accessed by decodePhrase */
393    saveInfo = info;
394    info->numberArguments = 0;
395    info->arguments = (char **) malloc(2 * sizeof(char *));
396    info->arguments[info->numberArguments++] = strdup("ampl");
397    info->arguments[info->numberArguments++] = strdup("cbc");
398    asl = ASL_alloc(ASL_read_f);
399    stub = getstub(&argv, &Oinfo);
400    if (!stub)
401        usage_ASL(&Oinfo, 1);
402    nl = jac0dim(stub, 0);
403    suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
404
405    /* set A_vals to get the constraints column-wise (malloc so can be freed) */
406    A_vals = (double *) malloc(nzc * sizeof(double));
407    if (!A_vals) {
408        printf("no memory\n");
409        return 1;
410    }
411    /* say we want primal solution */
412    want_xpi0 = 1;
413    /* for basis info */
414    info->columnStatus = (int *) malloc(n_var * sizeof(int));
415    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                    char ichar = sostype[i];
444                    assert (ichar == '1' || ichar == '2');
445                    info->sosType[i] = static_cast<char>(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            int 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        noNames_(false),
769        links_(0)
770{
771    problemName_ = "";
772    int status = 0;
773    if (!strcmp(fileName, "-") || !strcmp(fileName, "stdin")) {
774        // stdin
775    } else {
776        std::string name = fileName;
777        bool readable = fileCoinReadable(name);
778        if (!readable) {
779            std::cerr << "Unable to open file "
780                      << fileName << std::endl;
781            status = -1;
782        }
783    }
784    if (!status) {
785        gdb(nonLinear, fileName, info);
786    }
787}
788#ifdef JJF_ZERO
789static real
790qterm(ASL *asl, fint *colq, fint *rowq, real *delsq)
791{
792    double t, t1, *x, *x0, *xe;
793    fint *rq0, *rqe;
794
795    t = 0.;
796    x0 = x = X0;
797    xe = x + n_var;
798    rq0 = rowq;
799    while (x < xe) {
800        t1 = *x++;
801        rqe = rq0 + *++colq;
802        while (rowq < rqe)
803            t += t1 * x0[*rowq++]**delsq++;
804    }
805    return 0.5 * t;
806}
807#endif
808// stolen from IPopt with changes
809typedef struct {
810    double obj_sign_;
811    ASL_pfgh * asl_;
812    double * non_const_x_;
813    int * column_; // for jacobian
814    int * rowStart_;
815    double * gradient_;
816    double * constraintValues_;
817    int nz_h_full_; // number of nonzeros in hessian
818    int nerror_;
819    bool objval_called_with_current_x_;
820    bool conval_called_with_current_x_;
821    bool jacval_called_with_current_x_;
822} CbcAmplInfo;
823
824void
825CoinModel::gdb( int nonLinear, const char * fileName, const void * info)
826{
827    const ampl_info * amplInfo = (const ampl_info *) info;
828    ograd *og = NULL;
829    int i;
830    SufDesc *csd = NULL;
831    SufDesc *rsd = NULL;
832    /*bool *basis, *lower;*/
833    /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
834    //char tempBuffer[20];
835    double * objective = NULL;
836    double * columnLower = NULL;
837    double * columnUpper = NULL;
838    double * rowLower = NULL;
839    double * rowUpper = NULL;
840    int * columnStatus = NULL;
841    int * rowStatus = NULL;
842    int numberRows = -1;
843    int numberColumns = -1;
844    int numberElements = -1;
845    int numberBinary = -1;
846    int numberIntegers = -1;
847    int numberAllNonLinearBoth = 0;
848    int numberIntegerNonLinearBoth = 0;
849    int numberAllNonLinearConstraints = 0;
850    int numberIntegerNonLinearConstraints = 0;
851    int numberAllNonLinearObjective = 0;
852    int numberIntegerNonLinearObjective = 0;
853    double * primalSolution = NULL;
854    double direction = 1.0;
855    char * stub = strdup(fileName);
856    CoinPackedMatrix matrixByRow;
857    fint ** colqp = NULL;
858    int *z = NULL;
859    if (nonLinear == 0) {
860        // linear
861        asl = ASL_alloc(ASL_read_f);
862        nl = jac0dim(stub, 0);
863        free(stub);
864        suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
865
866        /* set A_vals to get the constraints column-wise (malloc so can be freed) */
867        A_vals = (double *) malloc(nzc * sizeof(double));
868        if (!A_vals) {
869            printf("no memory\n");
870            return ;
871        }
872        /* say we want primal solution */
873        want_xpi0 = 1;
874        /* for basis info */
875        columnStatus = (int *) malloc(n_var * sizeof(int));
876        rowStatus = (int *) malloc(n_con * sizeof(int));
877        csd = suf_iput("sstatus", ASL_Sufkind_var, columnStatus);
878        rsd = suf_iput("sstatus", ASL_Sufkind_con, rowStatus);
879        /* read linear model*/
880        f_read(nl, 0);
881        // see if any sos
882        if (true) {
883            char *sostype;
884            int nsosnz, *sosbeg, *sosind, * sospri;
885            double *sosref;
886            int nsos;
887            int i = ASL_suf_sos_explict_free;
888            int copri[2], **p_sospri;
889            copri[0] = 0;
890            copri[1] = 0;
891            p_sospri = &sospri;
892            nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
893                           &sosbeg, &sosind, &sosref);
894            if (nsos) {
895                abort();
896#ifdef JJF_ZERO
897                info->numberSos = nsos;
898                info->sosType = (char *) malloc(nsos);
899                info->sosPriority = (int *) malloc(nsos * sizeof(int));
900                info->sosStart = (int *) malloc((nsos + 1) * sizeof(int));
901                info->sosIndices = (int *) malloc(nsosnz * sizeof(int));
902                info->sosReference = (double *) malloc(nsosnz * sizeof(double));
903                sos_kludge(nsos, sosbeg, sosref, sosind);
904                for (int i = 0; i < nsos; i++) {
905                    int ichar = sostype[i];
906                    assert (ichar == '1' || ichar == '2');
907                    info->sosType[i] = ichar - '0';
908                }
909                memcpy(info->sosPriority, sospri, nsos*sizeof(int));
910                memcpy(info->sosStart, sosbeg, (nsos + 1)*sizeof(int));
911                memcpy(info->sosIndices, sosind, nsosnz*sizeof(int));
912                memcpy(info->sosReference, sosref, nsosnz*sizeof(double));
913#endif
914            }
915        }
916
917        /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
918        //Oinfo.uinfo = tempBuffer;
919        //if (getopts(argv, &Oinfo))
920        //return 1;
921        /* objective*/
922        objective = (double *) malloc(n_var * sizeof(double));
923        for (i = 0; i < n_var; i++)
924            objective[i] = 0.0;
925        if (n_obj) {
926            for (og = Ograd[0]; og; og = og->next)
927                objective[og->varno] = og->coef;
928        }
929        if (objtype[0])
930            direction = -1.0;
931        else
932            direction = 1.0;
933        objectiveOffset_ = objconst(0);
934        /* Column bounds*/
935        columnLower = (double *) malloc(n_var * sizeof(double));
936        columnUpper = (double *) malloc(n_var * sizeof(double));
937        for (i = 0; i < n_var; i++) {
938            columnLower[i] = LUv[2*i];
939            if (columnLower[i] <= negInfinity)
940                columnLower[i] = -COIN_DBL_MAX;
941            columnUpper[i] = LUv[2*i+1];
942            if (columnUpper[i] >= Infinity)
943                columnUpper[i] = COIN_DBL_MAX;
944        }
945        /* Row bounds*/
946        rowLower = (double *) malloc(n_con * sizeof(double));
947        rowUpper = (double *) malloc(n_con * sizeof(double));
948        for (i = 0; i < n_con; i++) {
949            rowLower[i] = LUrhs[2*i];
950            if (rowLower[i] <= negInfinity)
951                rowLower[i] = -COIN_DBL_MAX;
952            rowUpper[i] = LUrhs[2*i+1];
953            if (rowUpper[i] >= Infinity)
954                rowUpper[i] = COIN_DBL_MAX;
955        }
956        numberRows = n_con;
957        numberColumns = n_var;
958        numberElements = nzc;
959        numberBinary = nbv;
960        numberIntegers = niv;
961        /* put in primalSolution if exists */
962        if (X0) {
963            primalSolution = (double *) malloc(n_var * sizeof(double));
964            memcpy( primalSolution, X0, n_var*sizeof(double));
965        }
966        //double * dualSolution=NULL;
967        if (niv + nbv > 0)
968            mip_stuff(); // get any extra info
969        if ((!(niv + nbv) && (csd->kind & ASL_Sufkind_input))
970                || (rsd->kind & ASL_Sufkind_input)) {
971            /* convert status - need info on map */
972            static int map[] = {1, 3, 1, 1, 2, 1, 1};
973            stat_map(columnStatus, n_var, map, 6, "incoming columnStatus");
974            stat_map(rowStatus, n_con, map, 6, "incoming rowStatus");
975        } else {
976            /* all slack basis */
977            // leave status for output */
978#ifdef JJF_ZERO
979            free(rowStatus);
980            rowStatus = NULL;
981            free(columnStatus);
982            columnStatus = NULL;
983#endif
984        }
985        CoinPackedMatrix columnCopy(true, numberRows, numberColumns, numberElements,
986                                    A_vals, A_rownos, A_colstarts, NULL);
987        matrixByRow.reverseOrderedCopyOf(columnCopy);
988    } else if (nonLinear == 1) {
989        // quadratic
990        asl = ASL_alloc(ASL_read_fg);
991        nl = jac0dim(stub, (ftnlen) strlen(stub));
992        free(stub);
993        suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
994        /* read  model*/
995        X0 = (double*) malloc(n_var * sizeof(double));
996        CoinZeroN(X0, n_var);
997        qp_read(nl, 0);
998        assert (n_obj == 1);
999        int nz = 1 + n_con;
1000        colqp = (fint**) malloc(nz * (2 * sizeof(int*)
1001                                      + sizeof(double*)));
1002        fint ** rowqp = colqp + nz;
1003        double ** delsqp = (double **)(rowqp + nz);
1004        z = (int*) malloc(nz * sizeof(int));
1005        for (i = 0; i <= n_con; i++) {
1006            z[i] = nqpcheck(-i, rowqp + i, colqp + i, delsqp + i);
1007        }
1008        qp_opify();
1009        /* objective*/
1010        objective = (double *) malloc(n_var * sizeof(double));
1011        for (i = 0; i < n_var; i++)
1012            objective[i] = 0.0;
1013        if (n_obj) {
1014            for (og = Ograd[0]; og; og = og->next)
1015                objective[og->varno] = og->coef;
1016        }
1017        if (objtype[0])
1018            direction = -1.0;
1019        else
1020            direction = 1.0;
1021        objectiveOffset_ = objconst(0);
1022        /* Column bounds*/
1023        columnLower = (double *) malloc(n_var * sizeof(double));
1024        columnUpper = (double *) malloc(n_var * sizeof(double));
1025        for (i = 0; i < n_var; i++) {
1026            columnLower[i] = LUv[2*i];
1027            if (columnLower[i] <= negInfinity)
1028                columnLower[i] = -COIN_DBL_MAX;
1029            columnUpper[i] = LUv[2*i+1];
1030            if (columnUpper[i] >= Infinity)
1031                columnUpper[i] = COIN_DBL_MAX;
1032        }
1033        // Build by row from scratch
1034        //matrixByRow.reserve(n_var,nzc,true);
1035        // say row orderded
1036        matrixByRow.transpose();
1037        /* Row bounds*/
1038        rowLower = (double *) malloc(n_con * sizeof(double));
1039        rowUpper = (double *) malloc(n_con * sizeof(double));
1040        CoinBigIndex * rowStart = new CoinBigIndex [n_con+1];
1041        int * column = new int [nzc];
1042        double * element = new double [nzc];
1043        rowStart[0] = 0;
1044        numberElements = 0;
1045        for (i = 0; i < n_con; i++) {
1046            rowLower[i] = LUrhs[2*i];
1047            if (rowLower[i] <= negInfinity)
1048                rowLower[i] = -COIN_DBL_MAX;
1049            rowUpper[i] = LUrhs[2*i+1];
1050            if (rowUpper[i] >= Infinity)
1051                rowUpper[i] = COIN_DBL_MAX;
1052            for (cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
1053                column[numberElements] = cg->varno;
1054                element[numberElements++] = cg->coef;
1055            }
1056            rowStart[i+1] = numberElements;
1057        }
1058        assert (numberElements == nzc);
1059        matrixByRow.appendRows(n_con, rowStart, column, element);
1060        delete [] rowStart;
1061        delete [] column;
1062        delete [] element;
1063        numberRows = n_con;
1064        numberColumns = n_var;
1065        //numberElements=nzc;
1066        numberBinary = nbv;
1067        numberIntegers = niv;
1068        numberAllNonLinearBoth = nlvb;
1069        numberIntegerNonLinearBoth = nlvbi;
1070        numberAllNonLinearConstraints = nlvc;
1071        numberIntegerNonLinearConstraints = nlvci;
1072        numberAllNonLinearObjective = nlvo;
1073        numberIntegerNonLinearObjective = nlvoi;
1074        /* say we want primal solution */
1075        want_xpi0 = 1;
1076        //double * dualSolution=NULL;
1077        // save asl
1078        // Fix memory leak one day
1079        CbcAmplInfo * info = new CbcAmplInfo;
1080        //amplGamsData_ = info;
1081        info->asl_ = NULL; // as wrong form asl;
1082        info->nz_h_full_ = -1; // number of nonzeros in hessian
1083        info->objval_called_with_current_x_ = false;
1084        info->nerror_ = 0;
1085        info->obj_sign_ = direction;
1086        info->conval_called_with_current_x_ = false;
1087        info->non_const_x_ = NULL;
1088        info->jacval_called_with_current_x_ = false;
1089        info->rowStart_ = NULL;
1090        info->column_ = NULL;
1091        info->gradient_ = NULL;
1092        info->constraintValues_ = NULL;
1093    } else if (nonLinear == 2) {
1094        // General nonlinear!
1095        //ASL_pfgh* asl = (ASL_pfgh*)ASL_alloc(ASL_read_pfgh);
1096        asl = ASL_alloc(ASL_read_pfgh);
1097        nl = jac0dim(stub, (ftnlen) strlen(stub));
1098        free(stub);
1099        suf_declare(suftab, sizeof(suftab) / sizeof(SufDecl));
1100        /* read  model*/
1101        X0 = (double*) malloc(n_var * sizeof(double));
1102        CoinZeroN(X0, n_var);
1103        // code stolen from Ipopt
1104        int retcode = pfgh_read(nl, ASL_return_read_err | ASL_findgroups);
1105
1106        switch (retcode) {
1107        case ASL_readerr_none : {}
1108        break;
1109        case ASL_readerr_nofile : {
1110            printf( "Cannot open .nl file\n");
1111            exit(-1);
1112        }
1113        break;
1114        case ASL_readerr_nonlin : {
1115            assert(false); // this better not be an error!
1116            printf( "model involves nonlinearities (ed0read)\n");
1117            exit(-1);
1118        }
1119        break;
1120        case  ASL_readerr_argerr : {
1121            printf( "user-defined function with bad args\n");
1122            exit(-1);
1123        }
1124        break;
1125        case ASL_readerr_unavail : {
1126            printf( "user-defined function not available\n");
1127            exit(-1);
1128        }
1129        break;
1130        case ASL_readerr_corrupt : {
1131            printf( "corrupt .nl file\n");
1132            exit(-1);
1133        }
1134        break;
1135        case ASL_readerr_bug : {
1136            printf( "bug in .nl reader\n");
1137            exit(-1);
1138        }
1139        break;
1140        case ASL_readerr_CLP : {
1141            printf( "ASL error message: \"solver cannot handle CLP extensions\"\n");
1142            exit(-1);
1143        }
1144        break;
1145        default: {
1146            printf( "Unknown error in stub file read. retcode = %d\n", retcode);
1147            exit(-1);
1148        }
1149        break;
1150        }
1151
1152        // see "changes" in solvers directory of ampl code...
1153        hesset(1, 0, 1, 0, nlc);
1154
1155        assert (n_obj == 1);
1156        // find the nonzero structure for the hessian
1157        // parameters to sphsetup:
1158        int coeff_obj = 1; // coefficient of the objective fn ???
1159        int mult_supplied = 1; // multipliers will be supplied
1160        int uptri = 1; // only need the upper triangular part
1161        // save asl
1162        // Fix memory leak one day
1163        CbcAmplInfo * info = new CbcAmplInfo;
1164        moreInfo_ = (void *) info;
1165        //amplGamsData_ = info;
1166        info->asl_ = (ASL_pfgh *) asl;
1167        // This is not easy to get from ampl so save
1168        info->nz_h_full_ = sphsetup(-1, coeff_obj, mult_supplied, uptri);
1169        info->objval_called_with_current_x_ = false;
1170        info->nerror_ = 0;
1171        info->obj_sign_ = direction;
1172        info->conval_called_with_current_x_ = false;
1173        info->non_const_x_ = NULL;
1174        info->jacval_called_with_current_x_ = false;
1175        // Look at nonlinear
1176        if (nzc) {
1177            n_conjac[1] = nlc; // just nonlinear
1178            int * rowStart = new int [nlc+1];
1179            info->rowStart_ = rowStart;
1180            // See how many
1181            int  current_nz = 0;
1182            for (int i = 0; i < nlc; i++) {
1183                for (cgrad* cg = Cgrad[i]; cg; cg = cg->next) {
1184                    current_nz++;
1185                }
1186            }
1187            // setup the structure
1188            int * column = new int [current_nz];
1189            info->column_ = column;
1190            current_nz = 0;
1191            rowStart[0] = 0;
1192            for (int i = 0; i < nlc; i++) {
1193                for (cgrad* cg = Cgrad[i]; cg; cg = cg->next) {
1194                    cg->goff = current_nz;
1195                    //iRow[cg->goff] = i ;
1196                    //jCol[cg->goff] = cg->varno + 1;
1197                    column[cg->goff] = cg->varno ;
1198                    current_nz++;
1199                }
1200                rowStart[i+1] = current_nz;
1201            }
1202            info->gradient_ = new double [nzc];
1203            info->constraintValues_ = new double [nlc];
1204        }
1205        /* objective*/
1206        objective = (double *) malloc(n_var * sizeof(double));
1207        for (i = 0; i < n_var; i++)
1208            objective[i] = 0.0;
1209        if (n_obj) {
1210            for (og = Ograd[0]; og; og = og->next)
1211                objective[og->varno] = og->coef;
1212        }
1213        if (objtype[0])
1214            direction = -1.0;
1215        else
1216            direction = 1.0;
1217        objectiveOffset_ = objconst(0);
1218        /* Column bounds*/
1219        columnLower = (double *) malloc(n_var * sizeof(double));
1220        columnUpper = (double *) malloc(n_var * sizeof(double));
1221        for (i = 0; i < n_var; i++) {
1222            columnLower[i] = LUv[2*i];
1223            if (columnLower[i] <= negInfinity)
1224                columnLower[i] = -COIN_DBL_MAX;
1225            columnUpper[i] = LUv[2*i+1];
1226            if (columnUpper[i] >= Infinity)
1227                columnUpper[i] = COIN_DBL_MAX;
1228        }
1229        // Build by row from scratch
1230        //matrixByRow.reserve(n_var,nzc,true);
1231        // say row orderded
1232        matrixByRow.transpose();
1233        CoinBigIndex * rowStart = new CoinBigIndex [n_con+1];
1234        int * column = new int [nzc];
1235        double * element = new double [nzc];
1236        rowStart[0] = 0;
1237        numberElements = 0;
1238        /* Row bounds*/
1239        rowLower = (double *) malloc(n_con * sizeof(double));
1240        rowUpper = (double *) malloc(n_con * sizeof(double));
1241        for (i = 0; i < n_con; i++) {
1242            rowLower[i] = LUrhs[2*i];
1243            if (rowLower[i] <= negInfinity)
1244                rowLower[i] = -COIN_DBL_MAX;
1245            rowUpper[i] = LUrhs[2*i+1];
1246            if (rowUpper[i] >= Infinity)
1247                rowUpper[i] = COIN_DBL_MAX;
1248            for (cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
1249                column[numberElements] = cg->varno;
1250                double value = cg->coef;
1251                if (!value)
1252                    value = -1.2345e-29;
1253                element[numberElements++] = value;
1254            }
1255            rowStart[i+1] = numberElements;
1256        }
1257        assert (numberElements == nzc);
1258        matrixByRow.appendRows(n_con, rowStart, column, element);
1259        delete [] rowStart;
1260        delete [] column;
1261        delete [] element;
1262        numberRows = n_con;
1263        numberColumns = n_var;
1264        numberElements = nzc;
1265        numberBinary = nbv;
1266        numberIntegers = niv;
1267        numberAllNonLinearBoth = nlvb;
1268        numberIntegerNonLinearBoth = nlvbi;
1269        numberAllNonLinearConstraints = nlvc;
1270        numberIntegerNonLinearConstraints = nlvci;
1271        numberAllNonLinearObjective = nlvo;
1272        numberIntegerNonLinearObjective = nlvoi;
1273        /* say we want primal solution */
1274        want_xpi0 = 1;
1275        //double * dualSolution=NULL;
1276    } else {
1277        abort();
1278    }
1279    // set problem name
1280    problemName_ = "???";
1281
1282    // Build by row from scratch
1283    const double * element = matrixByRow.getElements();
1284    const int * column = matrixByRow.getIndices();
1285    const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
1286    const int * rowLength = matrixByRow.getVectorLengths();
1287    for (i = 0; i < numberRows; i++) {
1288        addRow(rowLength[i], column + rowStart[i],
1289               element + rowStart[i], rowLower[i], rowUpper[i]);
1290    }
1291    // Now do column part
1292    for (i = 0; i < numberColumns; i++) {
1293        setColumnBounds(i, columnLower[i], columnUpper[i]);
1294        setColumnObjective(i, objective[i]);
1295    }
1296    for ( i = numberColumns - numberBinary - numberIntegers;
1297            i < numberColumns; i++) {
1298        setColumnIsInteger(i, true);
1299    }
1300    // and non linear
1301    for (i = numberAllNonLinearBoth - numberIntegerNonLinearBoth;
1302            i < numberAllNonLinearBoth; i++) {
1303        setColumnIsInteger(i, true);
1304    }
1305    for (i = numberAllNonLinearConstraints - numberIntegerNonLinearConstraints;
1306            i < numberAllNonLinearConstraints; i++) {
1307        setColumnIsInteger(i, true);
1308    }
1309    for (i = numberAllNonLinearObjective - numberIntegerNonLinearObjective;
1310            i < numberAllNonLinearObjective; i++) {
1311        setColumnIsInteger(i, true);
1312    }
1313    free(columnLower);
1314    free(columnUpper);
1315    free(rowLower);
1316    free(rowUpper);
1317    free(objective);
1318    // do names
1319    int iRow;
1320    for (iRow = 0; iRow < numberRows_; iRow++) {
1321        char name[9];
1322        sprintf(name, "r%7.7d", iRow);
1323        setRowName(iRow, name);
1324    }
1325    int iColumn;
1326    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1327        char name[9];
1328        sprintf(name, "c%7.7d", iColumn);
1329        setColumnName(iColumn, name);
1330    }
1331    if (colqp) {
1332        // add in quadratic
1333        int nz = 1 + n_con;
1334        int nOdd = 0;
1335        fint ** rowqp = colqp + nz;
1336        double ** delsqp = (double **)(rowqp + nz);
1337        for (i = 0; i <= n_con; i++) {
1338            int nels = z[i];
1339            if (nels) {
1340                double * element = delsqp[i];
1341                int * start = (int *) colqp[i];
1342                int * row = (int *) rowqp[i];
1343                if (!element) {
1344                    // odd row - probably not quadratic
1345                    nOdd++;
1346                    continue;
1347                }
1348#ifdef JJF_ZERO
1349                printf("%d quadratic els\n", nels);
1350                for (int j = 0; j < n_var; j++) {
1351                    for (int k = start[j]; k < start[j+1]; k++)
1352                        printf("%d %d %g\n", j, row[k], element[k]);
1353                }
1354#endif
1355                if (i) {
1356                    int iRow = i - 1;
1357                    for (int j = 0; j < n_var; j++) {
1358                        for (int k = start[j]; k < start[j+1]; k++) {
1359                            int kColumn = row[k];
1360                            double value = element[k];
1361                            // ampl gives twice with assumed 0.5
1362                            if (kColumn < j)
1363                                continue;
1364                            else if (kColumn == j)
1365                                value *= 0.5;
1366                            const char * expr = getElementAsString(iRow, j);
1367                            double constant = 0.0;
1368                            bool linear;
1369                            if (expr && strcmp(expr, "Numeric")) {
1370                                linear = false;
1371                            } else {
1372                                constant = getElement(iRow, j);
1373                                linear = true;
1374                            }
1375                            char temp[1000];
1376                            char temp2[30];
1377                            if (value == 1.0)
1378                                sprintf(temp2, "c%7.7d", kColumn);
1379                            else
1380                                sprintf(temp2, "%g*c%7.7d", value, kColumn);
1381                            if (linear) {
1382                                if (!constant)
1383                                    strcpy(temp, temp2);
1384                                else if (value > 0.0)
1385                                    sprintf(temp, "%g+%s", constant, temp2);
1386                                else
1387                                    sprintf(temp, "%g%s", constant, temp2);
1388                            } else {
1389                                if (value > 0.0)
1390                                    sprintf(temp, "%s+%s", expr, temp2);
1391                                else
1392                                    sprintf(temp, "%s%s", expr, temp2);
1393                            }
1394                            assert (strlen(temp) < 1000);
1395                            setElement(iRow, j, temp);
1396                            if (amplInfo->logLevel > 1)
1397                                printf("el for row %d column c%7.7d is %s\n", iRow, j, temp);
1398                        }
1399                    }
1400                } else {
1401                    // objective
1402                    for (int j = 0; j < n_var; j++) {
1403                        for (int k = start[j]; k < start[j+1]; k++) {
1404                            int kColumn = row[k];
1405                            double value = element[k];
1406                            // ampl gives twice with assumed 0.5
1407                            if (kColumn < j)
1408                                continue;
1409                            else if (kColumn == j)
1410                                value *= 0.5;
1411                            const char * expr = getColumnObjectiveAsString(j);
1412                            double constant = 0.0;
1413                            bool linear;
1414                            if (expr && strcmp(expr, "Numeric")) {
1415                                linear = false;
1416                            } else {
1417                                constant = getColumnObjective(j);
1418                                linear = true;
1419                            }
1420                            char temp[1000];
1421                            char temp2[30];
1422                            if (value == 1.0)
1423                                sprintf(temp2, "c%7.7d", kColumn);
1424                            else
1425                                sprintf(temp2, "%g*c%7.7d", value, kColumn);
1426                            if (linear) {
1427                                if (!constant)
1428                                    strcpy(temp, temp2);
1429                                else if (value > 0.0)
1430                                    sprintf(temp, "%g+%s", constant, temp2);
1431                                else
1432                                    sprintf(temp, "%g%s", constant, temp2);
1433                            } else {
1434                                if (value > 0.0)
1435                                    sprintf(temp, "%s+%s", expr, temp2);
1436                                else
1437                                    sprintf(temp, "%s%s", expr, temp2);
1438                            }
1439                            assert (strlen(temp) < 1000);
1440                            setObjective(j, temp);
1441                            if (amplInfo->logLevel > 1)
1442                                printf("el for objective column c%7.7d is %s\n", j, temp);
1443                        }
1444                    }
1445                }
1446            }
1447        }
1448        if (nOdd) {
1449            printf("%d non-linear constraints could not be converted to quadratic\n", nOdd);
1450            exit(77);
1451        }
1452    }
1453    free(colqp);
1454    free(z);
1455    // see if any sos
1456    {
1457        char *sostype;
1458        int nsosnz, *sosbeg, *sosind, * sospri;
1459        double *sosref;
1460        int nsos;
1461        int i = ASL_suf_sos_explict_free;
1462        int copri[2], **p_sospri;
1463        copri[0] = 0;
1464        copri[1] = 0;
1465        p_sospri = &sospri;
1466        nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1467                       &sosbeg, &sosind, &sosref);
1468        if (nsos) {
1469            numberSOS_ = nsos;
1470            typeSOS_ = new int [numberSOS_];
1471            prioritySOS_ = new int [numberSOS_];
1472            startSOS_ = new int [numberSOS_+1];
1473            memberSOS_ = new int[nsosnz];
1474            referenceSOS_ = new double [nsosnz];
1475            sos_kludge(nsos, sosbeg, sosref, sosind);
1476            for (int i = 0; i < nsos; i++) {
1477                int ichar = sostype[i];
1478                assert (ichar == '1' || ichar == '2');
1479                typeSOS_[i] = ichar - '0';
1480            }
1481            memcpy(prioritySOS_, sospri, nsos*sizeof(int));
1482            memcpy(startSOS_, sosbeg, (nsos + 1)*sizeof(int));
1483            memcpy(memberSOS_, sosind, nsosnz*sizeof(int));
1484            memcpy(referenceSOS_, sosref, nsosnz*sizeof(double));
1485        }
1486    }
1487}
1488#else
1489#include "Cbc_ampl.h"
1490int
1491readAmpl(ampl_info * , int , char **, void ** )
1492{
1493    return 0;
1494}
1495void freeArrays1(ampl_info *)
1496{
1497}
1498void freeArrays2(ampl_info *)
1499{
1500}
1501void freeArgs(ampl_info * )
1502{
1503}
1504int ampl_obj_prec()
1505{
1506    return 0;
1507}
1508void writeAmpl(ampl_info * )
1509{
1510}
1511#endif
1512
Note: See TracBrowser for help on using the repository browser.