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

Last change on this file since 1370 was 1370, checked in by forrest, 11 years ago

add ids

  • Property svn:keywords set to Id
File size: 7.9 KB
Line 
1/* $Id: myPdco.cpp 1370 2009-06-04 09:37:13Z forrest $ */
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.