source: trunk/Clp/src/ClpLsqr.cpp @ 1304

Last change on this file since 1304 was 1171, checked in by forrest, 12 years ago

put back pdco in case wanted

File size: 11.3 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#ifdef COIN_DO_PDCO
4
5#include "ClpLsqr.hpp"
6#include "ClpPdco.hpp"
7
8
9void ClpLsqr::do_lsqr( CoinDenseVector<double> &b, 
10                double damp, double atol, double btol, double conlim, int itnlim, 
11                bool show, Info info, CoinDenseVector<double> &x , int *istop,
12                int *itn, Outfo *outfo, bool precon, CoinDenseVector<double> &Pr){
13
14 /**
15 Special version of LSQR for use with pdco.m.
16 It continues with a reduced atol if a pdco-specific test isn't
17 satisfied with the input atol.
18*/
19
20//     Initialize.
21  static char term_msg[8][80] = 
22  {
23    "The exact solution is x = 0",
24    "The residual Ax - b is small enough, given ATOL and BTOL",
25    "The least squares error is small enough, given ATOL",
26    "The estimated condition number has exceeded CONLIM",
27    "The residual Ax - b is small enough, given machine precision",
28    "The least squares error is small enough, given machine precision",
29    "The estimated condition number has exceeded machine precision",
30    "The iteration limit has been reached"
31  };
32
33  //  printf("***************** Entering LSQR *************\n");
34  assert (model_);
35 
36  char str1[100], str2[100], str3[100], str4[100], head1[100], head2[100];
37
38  int m = nrows_;  int n = ncols_;     // set  m,n from lsqr object
39 
40  *itn     = 0;            *istop   = 0;       int nstop  = 0;
41  double ctol   = 0;         
42  if (conlim > 0) ctol = 1/conlim;
43
44  double anorm  = 0;         double acond  = 0;
45  double dampsq = damp*damp; double ddnorm = 0;      double res2   = 0;
46  double xnorm  = 0;         double xxnorm = 0;      double z      = 0;
47  double cs2    = -1;        double sn2    = 0;
48
49  // Set up the first vectors u and v for the bidiagonalization.
50  // These satisfy  beta*u = b,  alfa*v = A'u.
51
52  CoinDenseVector<double> u(b);
53  CoinDenseVector<double> v(n, 0.0);
54  x.clear();
55  double alfa   = 0;      double beta = u.twoNorm();
56  if (beta > 0){
57    u = (1/beta) * u;
58    matVecMult( 2, v, u );
59    if (precon)
60      v = v*Pr;
61    alfa = v.twoNorm();
62  }
63  if (alfa > 0){
64    v.scale(1/alfa);
65  }
66  CoinDenseVector<double> w(v);
67
68  double arnorm = alfa * beta; 
69  if (arnorm == 0){
70     printf("  %s\n\n", term_msg[0]);
71     return;
72  }
73
74  double rhobar = alfa;  double phibar = beta;   double bnorm  = beta;
75  double rnorm  = beta;
76  sprintf(head1,"   Itn      x(1)      Function");
77  sprintf(head2," Compatible   LS      Norm A   Cond A");
78
79  if (show){
80    printf(" %s%s\n",head1, head2);
81    double test1  = 1;  double test2  = alfa / beta;
82    sprintf(str1,"%6d %12.5e %10.3e",   *itn, x[0], rnorm );
83    sprintf(str2,"  %8.1e  %8.1e",       test1, test2 );
84    printf("%s%s\n",str1, str2);
85  }
86
87  //----------------------------------------------------------------
88  // Main iteration loop.
89  //----------------------------------------------------------------
90  while (*itn < itnlim){
91    *itn += 1;
92    // Perform the next step of the bidiagonalization to obtain the
93    // next beta, u, alfa, v.  These satisfy the relations
94    // beta*u  =  a*v   -  alfa*u,
95    // alfa*v  =  A'*u  -  beta*v.
96
97    u.scale((-alfa));
98    if (precon){
99      CoinDenseVector<double> pv(v*Pr);
100      matVecMult( 1, u, pv);
101    }else{
102      matVecMult( 1, u, v);
103    }
104    beta = u.twoNorm();
105    if (beta > 0){
106      u.scale((1/beta));
107      anorm = sqrt(anorm*anorm + alfa*alfa + beta*beta +  damp*damp);
108      v.scale((-beta));
109      CoinDenseVector<double> vv(n);
110      vv.clear();
111      matVecMult( 2, vv, u );
112      if (precon)
113        vv = vv*Pr;
114      v = v + vv;
115      alfa  = v.twoNorm();
116      if (alfa > 0)
117        v.scale((1/alfa));
118    }
119
120    // Use a plane rotation to eliminate the damping parameter.
121    // This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
122
123    double rhobar1 = sqrt(rhobar*rhobar + damp*damp);
124    double cs1     = rhobar / rhobar1;
125    double sn1     = damp   / rhobar1;
126    double psi     = sn1 * phibar;
127    phibar       *= cs1;
128
129    // Use a plane rotation to eliminate the subdiagonal element (beta)
130    // of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
131
132    double rho     = sqrt(rhobar1*rhobar1 +  beta*beta);
133    double cs      =   rhobar1/ rho;
134    double sn      =   beta   / rho;
135    double theta   =   sn * alfa;
136    rhobar  = - cs * alfa;
137    double phi     =   cs * phibar;
138    phibar  =   sn * phibar;
139    double tau     =   sn * phi;
140
141    // Update x and w.
142
143    double t1      =   phi  /rho;
144    double t2      = - theta/rho;
145    //    dk           =   ((1/rho)*w);
146
147    double w_norm = w.twoNorm();
148    x       = x      +  t1*w;
149    w       = v      +  t2*w;
150    ddnorm  = ddnorm +  (w_norm/rho)*(w_norm/rho);
151    // if wantvar, var = var  +  dk.*dk; end
152
153    // Use a plane rotation on the right to eliminate the
154    // super-diagonal element (theta) of the upper-bidiagonal matrix.
155    // Then use the result to estimate  norm(x).
156
157    double delta   =   sn2 * rho;
158    double gambar  = - cs2 * rho;
159    double rhs     =   phi  -  delta * z;
160    double zbar    =   rhs / gambar;
161    xnorm   =   sqrt(xxnorm + zbar*zbar);
162    double gamma   =   sqrt(gambar*gambar + theta*theta);
163    cs2     =   gambar / gamma;
164    sn2     =   theta  / gamma;
165    z       =   rhs    / gamma;
166    xxnorm  =   xxnorm  +  z*z;
167
168    // Test for convergence.
169    // First, estimate the condition of the matrix  Abar,
170    // and the norms of  rbar  and  Abar'rbar.
171
172    acond   =   anorm * sqrt( ddnorm );
173    double res1    =   phibar*phibar;
174    double res2    =   res1  +  psi*psi;
175    rnorm   =   sqrt( res1 + res2 );
176    arnorm  =   alfa * fabs( tau );
177
178    // Now use these norms to estimate certain other quantities,
179    // some of which will be small near a solution.
180
181    double test1   =   rnorm / bnorm;
182    double test2   =   arnorm/( anorm * rnorm );
183    double test3   =       1 / acond;
184    t1      =   test1 / (1    +  anorm * xnorm / bnorm);
185    double rtol    =   btol  +  atol *  anorm * xnorm / bnorm;
186
187    // The following tests guard against extremely small values of
188    // atol, btol  or  ctol.  (The user may have set any or all of
189    // the parameters  atol, btol, conlim  to 0.)
190    // The effect is equivalent to the normal tests using
191    // atol = eps,  btol = eps,  conlim = 1/eps.
192
193    if (*itn >= itnlim)
194      *istop = 7; 
195    if (1 + test3  <= 1)
196      *istop = 6;
197    if (1 + test2  <= 1)
198      *istop = 5;
199    if (1 + t1     <= 1)
200      *istop = 4;
201
202    // Allow for tolerances set by the user.
203
204    if  (test3 <= ctol)
205      *istop = 3;
206    if  (test2 <= atol)
207      *istop = 2;
208    if  (test1 <= rtol)
209      *istop = 1;
210
211    //-------------------------------------------------------------------
212    // SPECIAL TEST THAT DEPENDS ON pdco.m.
213    // Aname in pdco   is  iw in lsqr.
214    // dy              is  x
215    // Other stuff     is in info.
216
217    // We allow for diagonal preconditioning in pdDDD3.
218    //-------------------------------------------------------------------
219    if (*istop > 0){
220      double r3new     = arnorm;
221      double r3ratio   = r3new / info.r3norm;
222      double atolold   = atol;
223      double atolnew   = atol;
224         
225      if (atol > info.atolmin){
226        if     (r3ratio <= 0.1){   // dy seems good
227          // Relax
228        }else if (r3ratio <= 0.5){ // Accept dy but make next one more accurate.
229          atolnew = atolnew * 0.1;
230        }else{                     // Recompute dy more accurately
231          if (show){
232            printf("\n                                ");
233            printf("                                \n");
234            printf(" %5.1f%7d%7.3f", log10(atolold), *itn, r3ratio);
235          }
236          atol    = atol * 0.1;
237          atolnew = atol;
238          *istop   = 0;
239        }
240   
241      outfo->atolold = atolold;
242      outfo->atolnew = atolnew;
243      outfo->r3ratio = r3ratio;
244    }
245     
246    //-------------------------------------------------------------------
247    // See if it is time to print something.
248    //-------------------------------------------------------------------
249    int prnt = 0;
250    if (n     <= 40       ) prnt = 1;
251    if (*itn   <= 10       ) prnt = 1;
252    if (*itn   >= itnlim-10) prnt = 1;
253    if (fmod(*itn,10) == 0  ) prnt = 1;
254    if (test3 <=  2*ctol  ) prnt = 1;
255    if (test2 <= 10*atol  ) prnt = 1;
256    if (test1 <= 10*rtol  ) prnt = 1;
257    if (*istop !=  0       ) prnt = 1;
258
259    if (prnt == 1){
260      if (show){
261        sprintf(str1, "   %6d %12.5e %10.3e",   *itn, x[0], rnorm );
262        sprintf(str2, "  %8.1e %8.1e",       test1, test2 );
263        sprintf(str3, " %8.1e %8.1e",        anorm, acond );
264        printf("%s%s%s\n", str1, str2, str3);
265      }
266    }
267    if (*istop > 0)
268      break;
269    }
270  }
271  // End of iteration loop.
272  // Print the stopping condition.
273
274  if (show){
275    printf("\n LSQR finished\n");
276  //    disp(msg(istop+1,:))
277  //    disp(' ')
278    printf("%s\n", term_msg[*istop]);
279    sprintf(str1, "istop  =%8d     itn    =%8d",      *istop, *itn    );
280    sprintf(str2, "anorm  =%8.1e   acond  =%8.1e",  anorm, acond  );
281    sprintf(str3, "rnorm  =%8.1e   arnorm =%8.1e",  rnorm, arnorm );
282    sprintf(str4, "bnorm  =%8.1e   xnorm  =%8.1e",  bnorm, xnorm  );
283    printf("%s %s\n", str1, str2);
284    printf("%s %s\n", str3, str4);
285  }
286}
287
288
289
290
291void ClpLsqr::matVecMult( int mode, CoinDenseVector<double> *x, CoinDenseVector<double> *y){
292  int n = model_->numberColumns();
293  int m = model_->numberRows();
294  CoinDenseVector<double> *temp = new CoinDenseVector<double>(n, 0.0);
295  double *t_elts = temp->getElements();
296  double *x_elts = x->getElements();
297  double *y_elts = y->getElements();
298  ClpPdco * pdcoModel = (ClpPdco *) model_;
299  if (mode == 1){
300    pdcoModel->matVecMult( 2, temp, y);
301    for (int k=0; k<n; k++)
302      x_elts[k] += (diag1_[k] * t_elts[k]);
303    for (int k=0; k<m; k++)
304      x_elts[n+k] += (diag2_ * y_elts[k]);
305  }else{
306    for (int k=0; k<n; k++)
307      t_elts[k] = diag1_[k] * y_elts[k];
308    pdcoModel->matVecMult( 1, x, temp);
309    for (int k=0; k<m; k++)
310      x_elts[k] += diag2_ * y_elts[n+k];
311  }
312  delete temp;
313  return;
314}
315
316void ClpLsqr::matVecMult( int mode, CoinDenseVector<double> &x, CoinDenseVector<double> &y){
317  matVecMult( mode, &x, &y);
318  return;
319}
320/* Default constructor */
321ClpLsqr::ClpLsqr() :
322  nrows_(0),
323  ncols_(0),
324  model_(NULL),
325  diag1_(NULL),
326  diag2_(0.0)
327{}
328
329/* Constructor for use with Pdco model (note modified for pdco!!!!) */
330ClpLsqr::ClpLsqr(ClpInterior *model) :
331  diag1_(NULL),
332  diag2_(0.0)
333{
334  model_ = model;
335  nrows_ = model->numberRows() + model->numberColumns();
336  ncols_ = model->numberRows();
337}
338/** Destructor */
339ClpLsqr::~ClpLsqr()
340{
341  // delete [] diag1_; no as we just borrowed it
342}
343bool 
344ClpLsqr::setParam(char *parmName, int parmValue){
345  std::cout << "Set lsqr integer parameter " << parmName << "to " << parmValue
346            << std::endl;
347  if ( strcmp(parmName, "nrows") == 0) {
348    nrows_ = parmValue;
349    return 1;
350  }else if ( strcmp(parmName, "ncols") == 0) {
351    ncols_ = parmValue;
352    return 1;
353  }
354  std::cout <<"Attempt to set unknown integer parameter name "<<parmName<<std::endl;
355  return 0;
356}
357ClpLsqr::ClpLsqr(const ClpLsqr &rhs) :
358  nrows_(rhs.nrows_),
359  ncols_(rhs.ncols_),
360  model_(rhs.model_),
361  diag2_(rhs.diag2_)
362{
363  diag1_ = ClpCopyOfArray(rhs.diag1_,nrows_);
364}
365// Assignment operator. This copies the data
366ClpLsqr & 
367ClpLsqr::operator=(const ClpLsqr & rhs)
368{
369  if (this != &rhs) {
370    delete [] diag1_;
371    diag1_ = ClpCopyOfArray(rhs.diag1_,nrows_);
372    nrows_ = rhs.nrows_;
373    ncols_ = rhs.ncols_;
374    model_ = rhs.model_;
375    diag2_ = rhs.diag2_;
376  }
377  return *this;
378}
379#endif
Note: See TracBrowser for help on using the repository browser.