source: stable/1.13/Clp/examples/myPdco.cpp @ 1898

Last change on this file since 1898 was 1559, checked in by stefan, 10 years ago

merge split branch into trunk

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