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

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

add ids

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