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

Last change on this file since 1665 was 1665, checked in by lou, 9 years ago

Add EPL license notice in src.

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