source: branches/pre/ClpLsqr.cpp @ 1926

Last change on this file since 1926 was 219, checked in by forrest, 17 years ago

pdco seems to work

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