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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

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