source: trunk/Clp/examples/myPdco.cpp @ 1362

Last change on this file since 1362 was 1276, checked in by ladanyi, 12 years ago

changed max/min to CoinMax/CoinMin?

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