source: trunk/Clp/examples/testQP.cpp @ 778

Last change on this file since 778 was 778, checked in by andreasw, 13 years ago

corrected Externals to exclude MSDevStudio; made Clp examples to work

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