source: trunk/Clp/examples/testQP.cpp

Last change on this file was 2278, checked in by forrest, 4 months ago

COIN_BIG_INDEX 2 changes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 7.6 KB
Line 
1/* $Id: testQP.cpp 2278 2017-10-02 09:51:14Z forrest $ */
2// Copyright (C) 2005, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinMpsIO.hpp"
7#include "ClpInterior.hpp"
8#include "ClpSimplex.hpp"
9#include "ClpCholeskyBase.hpp"
10#include "ClpQuadraticObjective.hpp"
11#include <cassert>
12int main(int argc, const char *argv[])
13{
14     /* Read quadratic model in two stages to test loadQuadraticObjective.
15
16        And is also possible to just read into ClpSimplex/Interior which sets it all up in one go.
17        But this is only if it is in QUADOBJ format.
18
19        If no arguments does share2qp using ClpInterior (also creates quad.mps which is in QUADOBJ format)
20        If one argument uses simplex e.g. testit quad.mps
21        If > one uses barrier via ClpSimplex input and then ClpInterior borrow
22     */
23     if (argc < 2) {
24          CoinMpsIO  m;
25#if defined(SAMPLEDIR)
26          int status = m.readMps(SAMPLEDIR "/share2qp", "mps");
27          if (status) {
28               printf("errors on input\n");
29               exit(77);
30          }
31#else
32          fprintf(stderr, "Do not know where to find sample MPS files.\n");
33          exit(1);
34#endif
35          ClpInterior model;
36          model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
37                            m.getObjCoefficients(),
38                            m.getRowLower(), m.getRowUpper());
39          // get quadratic part
40          CoinBigIndex * start = NULL;
41          int * column = NULL;
42          double * element = NULL;
43          m.readQuadraticMps(NULL, start, column, element, 2);
44          int j;
45          for (j = 0; j < 79; j++) {
46               if (start[j] < start[j+1]) {
47                    CoinBigIndex i;
48                    printf("Column %d ", j);
49                    for (i = start[j]; i < start[j+1]; i++) {
50                         printf("( %d, %g) ", column[i], element[i]);
51                    }
52                    printf("\n");
53               }
54          }
55          model.loadQuadraticObjective(model.numberColumns(), start, column, element);
56          // share2qp is in old style qp format - convert to new so other options can use
57          model.writeMps("quad.mps");
58          ClpCholeskyBase * cholesky = new ClpCholeskyBase();
59          cholesky->setKKT(true);
60          model.setCholesky(cholesky);
61          model.primalDual();
62          double *primal;
63          double *dual;
64          primal = model.primalColumnSolution();
65          dual = model.dualRowSolution();
66          int i;
67          int numberColumns = model.numberColumns();
68          int numberRows = model.numberRows();
69          for (i = 0; i < numberColumns; i++) {
70               if (fabs(primal[i]) > 1.0e-8)
71                    printf("%d primal %g\n", i, primal[i]);
72          }
73          for (i = 0; i < numberRows; i++) {
74               if (fabs(dual[i]) > 1.0e-8)
75                    printf("%d dual %g\n", i, dual[i]);
76          }
77     } else {
78          // Could read into ClpInterior
79          ClpSimplex model;
80          if (model.readMps(argv[1])) {
81               printf("errors on input\n");
82               exit(77);
83          }
84          model.writeMps("quad");
85          if (argc < 3) {
86               // simplex - just primal as dual does not work
87               // also I need to fix scaling of duals on output
88               // (Was okay in first place - can't mix and match scaling techniques)
89               // model.scaling(0);
90               model.primal();
91          } else {
92               // barrier
93               ClpInterior barrier;
94               barrier.borrowModel(model);
95               ClpCholeskyBase * cholesky = new ClpCholeskyBase();
96               cholesky->setKKT(true);
97               barrier.setCholesky(cholesky);
98               barrier.primalDual();
99               barrier.returnModel(model);
100          }
101          // Just check if share2qp (quad.mps here)
102          // this is because I am not checking if variables at ub
103          if (model.numberColumns() == 79) {
104               double *primal;
105               double *dual;
106               primal = model.primalColumnSolution();
107               dual = model.dualRowSolution();
108               // Check duals by hand
109               const ClpQuadraticObjective * quadraticObj =
110                    (dynamic_cast<const ClpQuadraticObjective*>(model.objectiveAsObject()));
111               assert(quadraticObj);
112               CoinPackedMatrix * quad = quadraticObj->quadraticObjective();
113               const int * columnQuadratic = quad->getIndices();
114               const CoinBigIndex * columnQuadraticStart = quad->getVectorStarts();
115               const int * columnQuadraticLength = quad->getVectorLengths();
116               const double * quadraticElement = quad->getElements();
117               int numberColumns = model.numberColumns();
118               int numberRows = model.numberRows();
119               double * gradient = new double [numberColumns];
120               // move linear objective
121               memcpy(gradient, quadraticObj->linearObjective(), numberColumns * sizeof(double));
122               int iColumn;
123               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
124                    double valueI = primal[iColumn];
125                    CoinBigIndex j;
126                    for (j = columnQuadraticStart[iColumn];
127                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
128                         int jColumn = columnQuadratic[j];
129                         double valueJ = primal[jColumn];
130                         double elementValue = quadraticElement[j];
131                         if (iColumn != jColumn) {
132                              double gradientI = valueJ * elementValue;
133                              double gradientJ = valueI * elementValue;
134                              gradient[iColumn] += gradientI;
135                              gradient[jColumn] += gradientJ;
136                         } else {
137                              double gradientI = valueI * elementValue;
138                              gradient[iColumn] += gradientI;
139                         }
140                    }
141                    if (fabs(primal[iColumn]) > 1.0e-8)
142                         printf("%d primal %g\n", iColumn, primal[iColumn]);
143               }
144               for (int i = 0; i < numberRows; i++) {
145                    if (fabs(dual[i]) > 1.0e-8)
146                         printf("%d dual %g\n", i, dual[i]);
147               }
148               // Now use duals to get reduced costs
149               // Can't use this as will try and use scaling
150               // model.transposeTimes(-1.0,dual,gradient);
151               // So ...
152               CoinPackedMatrix * matrix = model.matrix();
153               const int * row = matrix->getIndices();
154               const CoinBigIndex * columnStart = matrix->getVectorStarts();
155               const int * columnLength = matrix->getVectorLengths();
156               const double * element = matrix->getElements();
157               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
158                    double dj = gradient[iColumn];
159                    CoinBigIndex j;
160                    for (j = columnStart[iColumn];
161                              j < columnStart[iColumn] + columnLength[iColumn]; j++) {
162                         int jRow = row[j];
163                         dj -= element[j] * dual[jRow];
164                    }
165                    if (model.getColumnStatus(iColumn) == ClpSimplex::basic) {
166                         assert(fabs(dj) < 1.0e-5);
167                    } else {
168                         assert(dj > -1.0e-5);
169                    }
170               }
171               delete [] gradient;
172          }
173     }
174     return 0;
175}
Note: See TracBrowser for help on using the repository browser.