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

Last change on this file since 2385 was 2385, checked in by unxusr, 4 months ago

formatting

  • Property svn:keywords set to Id
File size: 11.2 KB
Line 
1/* $Id: ClpLsqr.cpp 2385 2019-01-06 19:43:06Z unxusr $ */
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
9void ClpLsqr::do_lsqr(CoinDenseVector< double > &b,
10  double damp, double atol, double btol, double conlim, int itnlim,
11  bool show, Info info, CoinDenseVector< double > &x, int *istop,
12  int *itn, Outfo *outfo, bool precon, CoinDenseVector< double > &Pr)
13{
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    "The exact solution is x = 0",
24    "The residual Ax - b is small enough, given ATOL and BTOL",
25    "The least squares error is small enough, given ATOL",
26    "The estimated condition number has exceeded CONLIM",
27    "The residual Ax - b is small enough, given machine precision",
28    "The least squares error is small enough, given machine precision",
29    "The estimated condition number has exceeded machine precision",
30    "The iteration limit has been reached"
31  };
32
33  //  printf("***************** Entering LSQR *************\n");
34  assert(model_);
35
36  char str1[100], str2[100], str3[100], str4[100], head1[100], head2[100];
37
38  int n = ncols_; // set  m,n from lsqr object
39
40  *itn = 0;
41  *istop = 0;
42  double ctol = 0;
43  if (conlim > 0)
44    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)
261        prnt = 1;
262      if (*itn <= 10)
263        prnt = 1;
264      if (*itn >= itnlim - 10)
265        prnt = 1;
266      if (*itn % 10 == 0)
267        prnt = 1;
268      if (test3 <= 2 * ctol)
269        prnt = 1;
270      if (test2 <= 10 * atol)
271        prnt = 1;
272      if (test1 <= 10 * rtol)
273        prnt = 1;
274      if (*istop != 0)
275        prnt = 1;
276
277      if (prnt == 1) {
278        if (show) {
279          sprintf(str1, "   %6d %12.5e %10.3e", *itn, x[0], rnorm);
280          sprintf(str2, "  %8.1e %8.1e", test1, test2);
281          sprintf(str3, " %8.1e %8.1e", anorm, acond);
282          printf("%s%s%s\n", str1, str2, str3);
283        }
284      }
285      if (*istop > 0)
286        break;
287    }
288  }
289  // End of iteration loop.
290  // Print the stopping condition.
291
292  if (show) {
293    printf("\n LSQR finished\n");
294    //    disp(msg(istop+1,:))
295    //    disp(' ')
296    printf("%s\n", term_msg[*istop]);
297    sprintf(str1, "istop  =%8d     itn    =%8d", *istop, *itn);
298    sprintf(str2, "anorm  =%8.1e   acond  =%8.1e", anorm, acond);
299    sprintf(str3, "rnorm  =%8.1e   arnorm =%8.1e", rnorm, arnorm);
300    sprintf(str4, "bnorm  =%8.1e   xnorm  =%8.1e", bnorm, xnorm);
301    printf("%s %s\n", str1, str2);
302    printf("%s %s\n", str3, str4);
303  }
304}
305
306void ClpLsqr::matVecMult(int mode, CoinDenseVector< double > *x, CoinDenseVector< double > *y)
307{
308  int n = model_->numberColumns();
309  int m = model_->numberRows();
310  CoinDenseVector< double > *temp = new CoinDenseVector< double >(n, 0.0);
311  double *t_elts = temp->getElements();
312  double *x_elts = x->getElements();
313  double *y_elts = y->getElements();
314  ClpPdco *pdcoModel = (ClpPdco *)model_;
315  if (mode == 1) {
316    pdcoModel->matVecMult(2, temp, y);
317    for (int k = 0; k < n; k++)
318      x_elts[k] += (diag1_[k] * t_elts[k]);
319    for (int k = 0; k < m; k++)
320      x_elts[n + k] += (diag2_ * y_elts[k]);
321  } else {
322    for (int k = 0; k < n; k++)
323      t_elts[k] = diag1_[k] * y_elts[k];
324    pdcoModel->matVecMult(1, x, temp);
325    for (int k = 0; k < m; k++)
326      x_elts[k] += diag2_ * y_elts[n + k];
327  }
328  delete temp;
329  return;
330}
331
332void ClpLsqr::matVecMult(int mode, CoinDenseVector< double > &x, CoinDenseVector< double > &y)
333{
334  matVecMult(mode, &x, &y);
335  return;
336}
337/* Default constructor */
338ClpLsqr::ClpLsqr()
339  : nrows_(0)
340  , ncols_(0)
341  , model_(NULL)
342  , diag1_(NULL)
343  , diag2_(0.0)
344{
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 ClpLsqr::setParam(char *parmName, int parmValue)
362{
363  std::cout << "Set lsqr integer parameter " << parmName << "to " << parmValue
364            << std::endl;
365  if (strcmp(parmName, "nrows") == 0) {
366    nrows_ = parmValue;
367    return 1;
368  } else if (strcmp(parmName, "ncols") == 0) {
369    ncols_ = parmValue;
370    return 1;
371  }
372  std::cout << "Attempt to set unknown integer parameter name " << parmName << std::endl;
373  return 0;
374}
375ClpLsqr::ClpLsqr(const ClpLsqr &rhs)
376  : nrows_(rhs.nrows_)
377  , ncols_(rhs.ncols_)
378  , model_(rhs.model_)
379  , diag2_(rhs.diag2_)
380{
381  diag1_ = ClpCopyOfArray(rhs.diag1_, nrows_);
382}
383// Assignment operator. This copies the data
384ClpLsqr &
385ClpLsqr::operator=(const ClpLsqr &rhs)
386{
387  if (this != &rhs) {
388    delete[] diag1_;
389    diag1_ = ClpCopyOfArray(rhs.diag1_, nrows_);
390    nrows_ = rhs.nrows_;
391    ncols_ = rhs.ncols_;
392    model_ = rhs.model_;
393    diag2_ = rhs.diag2_;
394  }
395  return *this;
396}
397
398/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
399*/
Note: See TracBrowser for help on using the repository browser.