source: trunk/Clp/examples/myPdco.cpp

Last change on this file was 1662, checked in by lou, 8 years ago

Add EPL license notice in examples.

  • Property svn:keywords set to Id
File size: 9.0 KB
Line 
1/* $Id: myPdco.cpp 1662 2011-01-04 17:52:40Z forrest $ */
2// Copyright (C) 2003, 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
7#include <cstdio>
8
9#include "CoinPragma.hpp"
10#include "CoinHelperFunctions.hpp"
11
12#include "ClpInterior.hpp"
13#include "myPdco.hpp"
14#include "ClpDummyMatrix.hpp"
15#include "ClpMessage.hpp"
16
17//#############################################################################
18// Constructors / Destructor / Assignment
19//#############################################################################
20
21//-------------------------------------------------------------------
22// Default Constructor
23//-------------------------------------------------------------------
24myPdco::myPdco()
25     : ClpPdcoBase(),
26       rowIndex_(NULL),
27       numlinks_(0),
28       numnodes_(0)
29{
30     setType(11);
31}
32
33// Constructor from stuff
34myPdco::myPdco(double d1, double d2,
35               int numnodes, int numlinks)
36     : ClpPdcoBase(),
37       rowIndex_(NULL),
38       numlinks_(numlinks),
39       numnodes_(numnodes)
40{
41     d1_ = d1;
42     d2_ = d2;
43     setType(11);
44}
45//-------------------------------------------------------------------
46// Copy constructor
47//-------------------------------------------------------------------
48myPdco::myPdco(const myPdco & rhs)
49     : ClpPdcoBase(rhs),
50       numlinks_(rhs.numlinks_),
51       numnodes_(rhs.numnodes_)
52{
53     rowIndex_ = ClpCopyOfArray(rhs.rowIndex_, 2 * (numlinks_ + 2 * numnodes_));
54}
55
56//-------------------------------------------------------------------
57// Destructor
58//-------------------------------------------------------------------
59myPdco::~myPdco()
60{
61     delete [] rowIndex_;
62}
63
64//----------------------------------------------------------------
65// Assignment operator
66//-------------------------------------------------------------------
67myPdco &
68myPdco::operator= (const myPdco & rhs)
69{
70     if (this != &rhs) {
71          ClpPdcoBase::operator= (rhs);
72          numlinks_ = rhs.numlinks_;
73          numnodes_ = rhs.numnodes_;
74          rowIndex_ = ClpCopyOfArray(rhs.rowIndex_, 2 * (numlinks_ + 2 * numnodes_));
75     }
76     return *this;
77}
78//-------------------------------------------------------------------
79// Clone
80//-------------------------------------------------------------------
81ClpPdcoBase * myPdco::clone() const
82{
83     return new myPdco(*this);
84}
85
86void myPdco::matVecMult(ClpInterior * model,  int mode, double* x_elts, double* y_elts) const
87{
88     int nrow = model->numberRows();
89     if (mode == 1) {
90          double y_sum = 0.0;
91          for (int k = 0; k < numlinks_; k++) {
92               y_sum += y_elts[k];
93               int i1 = rowIndex_[2*k];
94               assert(i1 >= 0);
95               x_elts[i1] += y_elts[k];
96               int i2 = rowIndex_[2*k+1];
97               assert(i2 >= 0);
98               x_elts[i2] -= y_elts[k];
99          }
100          double y_suma = 0.0;
101          double y_sumb = 0.0;
102          for (int k = 0; k < numnodes_; k++) {
103               y_suma -= y_elts[numlinks_ + k];
104               y_sumb += y_elts[numlinks_ + numnodes_ + k];
105               x_elts[k] += y_elts[numlinks_ + k];
106               x_elts[k] -= y_elts[numlinks_ + numnodes_ + k];
107          }
108          x_elts[nrow-3] += y_suma;
109          x_elts[nrow-2] += y_sumb;
110          x_elts[nrow-1] += (y_sum - y_suma + y_sumb);
111     } else {
112          for (int k = 0; k < numlinks_; k++) {
113               x_elts[k] += y_elts[nrow-1];
114               int i1 = rowIndex_[2*k];
115               x_elts[k] += y_elts[i1];
116               int i2 = rowIndex_[2*k+1];
117               x_elts[k] -= y_elts[i2];
118          }
119          for (int k = 0; k < numnodes_; k++) {
120               x_elts[numlinks_ + k] += (y_elts[k] - y_elts[nrow-3] + y_elts[nrow-1]);
121               x_elts[numlinks_ + numnodes_ + k] +=
122                    (y_elts[nrow-2] - y_elts[k] + y_elts[nrow-1]);
123          }
124     }
125     return;
126}
127
128void myPdco::matPrecon(ClpInterior * model, double delta, double* x_elts, double* y_elts) const
129{
130     double y_sum = 0.0;
131     int ncol = model->numberColumns();
132     int nrow = model->numberRows();
133     double *ysq = new double[ncol];
134     for (int k = 0; k < nrow; k++)
135          x_elts[k] = 0.0;
136     for (int k = 0; k < ncol; k++)
137          ysq[k] = y_elts[k] * y_elts[k];
138
139     for (int k = 0; k < numlinks_; k++) {
140          y_sum += ysq[k];
141          int i1 = rowIndex_[2*k];
142          x_elts[i1] += ysq[k];
143          int i2 = rowIndex_[2*k+1];
144          x_elts[i2] += ysq[k];
145     }
146     double y_suma = 0.0;
147     double y_sumb = 0.0;
148     for (int k = 0; k < numnodes_; k++) {
149          y_suma += ysq[numlinks_ + k];
150          y_sumb += ysq[numlinks_ + numnodes_ + k];
151          x_elts[k] += ysq[numlinks_ + k];
152          x_elts[k] += ysq[numlinks_ + numnodes_ + k];
153     }
154     x_elts[nrow-3] += y_suma;
155     x_elts[nrow-2] += y_sumb;
156     x_elts[nrow-1] += (y_sum + y_suma + y_sumb);
157     delete [] ysq;
158     double delsq = delta * delta;
159     for (int k = 0; k < nrow; k++)
160          x_elts[k] = 1.0 / (sqrt(x_elts[k] + delsq));
161     return;
162}
163double myPdco::getObj(ClpInterior * model, CoinDenseVector<double> &x) const
164{
165     double obj = 0;
166     double *x_elts = x.getElements();
167     int ncol = model->numberColumns();
168     for (int k = 0; k < ncol; k++)
169          obj += x_elts[k] * log(x_elts[k]);
170     return obj;
171}
172
173void myPdco::getGrad(ClpInterior * model, CoinDenseVector<double> &x, CoinDenseVector<double> &g) const
174{
175     double *x_elts = x.getElements();
176     double *g_elts = g.getElements();
177     int ncol = model->numberColumns();
178     for (int k = 0; k < ncol; k++)
179          g_elts[k] = 1.0 + log(x_elts[k]);
180     return;
181}
182
183void myPdco::getHessian(ClpInterior * model, CoinDenseVector<double> &x, CoinDenseVector<double> &H) const
184{
185     double *x_elts = x.getElements();
186     double *H_elts = H.getElements();
187     int ncol = model->numberColumns();
188     for (int k = 0; k < ncol; k++)
189          H_elts[k] = 1.0 / x_elts[k];
190     return;
191}
192myPdco::myPdco(ClpInterior & model, FILE * fpData, FILE * fpParam)
193{
194     int nrow;
195     int ncol;
196     int numelts;
197     double  *rowUpper;
198     double  *rowLower;
199     double  *colUpper;
200     double  *colLower;
201     double  *rhs;
202     double  *x;
203     double  *y;
204     double  *dj;
205     int ipair[2], igparm[4], maxrows, maxlinks;
206     // Read max array sizes and allocate them
207     fscanf(fpParam, "%d %d", &maxrows, &maxlinks);
208     int *ir = new int[2*maxlinks + 5];
209     /********************** alpha parameter hrdwired here ***********/
210     double alpha = 0.9;
211
212     int kct = 0;
213     int imax = 0;
214     int imin = 0x7fffffff;
215     int *ifrom = &ipair[0];
216     int *ito = &ipair[1];
217     int nonzpt = 0;
218     numlinks_ = 0;
219     while (fscanf(fpData, "%d %d", ifrom, ito) && kct++ < maxlinks) {
220          //  while(fread(ipair, 4,2, fpData) && kct++ < maxlinks){
221          if ((*ifrom) < 0) {
222               printf("Bad link  %d  %d\n", *ifrom, *ito);
223               continue;
224          }
225          ipair[0]--;
226          ipair[1]--;
227          //assert(*ifrom>=0&&*ifrom<maxrows);
228          //assert(*ito>=0&&*ifrom<maxrows);
229          if (*ifrom < 0 || *ifrom >= maxrows || *ito < 0 || *ito >= maxrows) {
230               printf("bad link %d %d\n", *ifrom, *ito);
231               continue;
232          }
233          numlinks_++;
234          ir[nonzpt++] = *ifrom;
235          ir[nonzpt++] = *ito;
236          imax = CoinMax(imax, *ifrom);
237          imax = CoinMax(imax, *ito);
238          imin = CoinMin(imin, *ifrom);
239          imin = CoinMin(imin, *ito);
240     }
241     fclose(fpData);
242     fclose(fpParam);
243     printf("imax %d  imin %d\n", imax, imin);
244     // Set model size
245     numnodes_ = imax + 1;
246     nrow = numnodes_ + 3;
247     ncol = numlinks_ + 2 * numnodes_;
248     numelts = 3 * ncol;
249
250     rowIndex_ = ir;
251
252     d1_ = 1.0e-3;
253     d2_ = 1.0e-3;
254
255     double* rhs_def = new double[nrow];
256     for (int k = 0; k < nrow; k++)
257          rhs_def[k] = 0.0;
258     rhs_def[nrow-3] = alpha - 1.0;
259     rhs_def[nrow-2] = 1.0 - alpha;
260     rhs_def[nrow-1] = 1.0;
261     // rhs_ etc should not be public
262     rhs = rhs_def;
263     rowUpper = rhs_def;
264     rowLower = rhs_def;
265     double *x_def = new double[ncol];
266     double *U_def = new double[ncol];
267     double *L_def = new double[ncol];
268     for (int k = 0; k < ncol; k++) {
269          //    x_def[k] = 10.0/ncol;
270          x_def[k] = 1.0 / ncol;
271          U_def[k] = 1e20;;
272          L_def[k] = 0.0;
273     }
274     x = x_def;
275     colUpper = U_def;
276     colLower = L_def;
277     // We have enough to create a model
278     ClpDummyMatrix dummy(ncol, nrow, numelts);
279     model.loadProblem(dummy,
280                       colLower, colUpper, NULL, rowLower, rowUpper);
281     double *y_def = new double[nrow];
282     for (int k = 0; k < nrow; k++)
283          y_def[k] = 0.0;
284     y = y_def;
285     double *dj_def = new double[ncol];
286     for (int k = 0; k < ncol; k++)
287          dj_def[k] = 1.0;
288     dj = dj_def;
289     // delete arrays
290     delete [] U_def;
291     delete [] L_def;
292     // Should be sets
293     model.rhs_ = rhs;
294     model.x_ = x;
295     model.y_ = y;
296     model.dj_ = dj;
297     model.xsize_ = 50 / ncol;
298     model.xsize_ = CoinMin(model.xsize_, 1.0);
299     model.zsize_ = 1;
300}
Note: See TracBrowser for help on using the repository browser.