source: stable/1.15/Clp/src/ClpLsqr.cpp @ 1949

Last change on this file since 1949 was 1941, checked in by stefan, 7 years ago

sync with trunk rev 1940

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