source: trunk/Clp/src/Clp_ampl.cpp @ 2030

Last change on this file since 2030 was 2030, checked in by forrest, 6 years ago

fix some ampl stuff, add ClpSolver? and a few fixes

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