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
399