1 | /* $Id: ClpLsqr.hpp 2385 2019-01-06 19:43:06Z stefan $ */ |
---|
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 | #ifndef ClpLsqr_H_ |
---|
7 | #define ClpLsqr_H_ |
---|
8 | |
---|
9 | #include "CoinDenseVector.hpp" |
---|
10 | |
---|
11 | #include "ClpInterior.hpp" |
---|
12 | |
---|
13 | /** |
---|
14 | This class implements LSQR |
---|
15 | |
---|
16 | @verbatim |
---|
17 | LSQR solves Ax = b or min ||b - Ax||_2 if damp = 0, |
---|
18 | or min || (b) - ( A )x || otherwise. |
---|
19 | || (0) (damp I) ||2 |
---|
20 | A is an m by n matrix defined by user provided routines |
---|
21 | matVecMult(mode, y, x) |
---|
22 | which performs the matrix-vector operations where y and x |
---|
23 | are references or pointers to CoinDenseVector objects. |
---|
24 | If mode = 1, matVecMult must return y = Ax without altering x. |
---|
25 | If mode = 2, matVecMult must return y = A'x without altering x. |
---|
26 | |
---|
27 | ----------------------------------------------------------------------- |
---|
28 | LSQR uses an iterative (conjugate-gradient-like) method. |
---|
29 | For further information, see |
---|
30 | 1. C. C. Paige and M. A. Saunders (1982a). |
---|
31 | LSQR: An algorithm for sparse linear equations and sparse least squares, |
---|
32 | ACM TOMS 8(1), 43-71. |
---|
33 | 2. C. C. Paige and M. A. Saunders (1982b). |
---|
34 | Algorithm 583. LSQR: Sparse linear equations and least squares problems, |
---|
35 | ACM TOMS 8(2), 195-209. |
---|
36 | 3. M. A. Saunders (1995). Solution of sparse rectangular systems using |
---|
37 | LSQR and CRAIG, BIT 35, 588-604. |
---|
38 | |
---|
39 | Input parameters: |
---|
40 | |
---|
41 | atol, btol are stopping tolerances. If both are 1.0e-9 (say), |
---|
42 | the final residual norm should be accurate to about 9 digits. |
---|
43 | (The final x will usually have fewer correct digits, |
---|
44 | depending on cond(A) and the size of damp.) |
---|
45 | conlim is also a stopping tolerance. lsqr terminates if an estimate |
---|
46 | of cond(A) exceeds conlim. For compatible systems Ax = b, |
---|
47 | conlim could be as large as 1.0e+12 (say). For least-squares |
---|
48 | problems, conlim should be less than 1.0e+8. |
---|
49 | Maximum precision can be obtained by setting |
---|
50 | atol = btol = conlim = zero, but the number of iterations |
---|
51 | may then be excessive. |
---|
52 | itnlim is an explicit limit on iterations (for safety). |
---|
53 | show = 1 gives an iteration log, |
---|
54 | show = 0 suppresses output. |
---|
55 | info is a structure special to pdco.m, used to test if |
---|
56 | was small enough, and continuing if necessary with smaller atol. |
---|
57 | |
---|
58 | |
---|
59 | Output parameters: |
---|
60 | x is the final solution. |
---|
61 | *istop gives the reason for termination. |
---|
62 | *istop = 1 means x is an approximate solution to Ax = b. |
---|
63 | = 2 means x approximately solves the least-squares problem. |
---|
64 | rnorm = norm(r) if damp = 0, where r = b - Ax, |
---|
65 | = sqrt( norm(r)**2 + damp**2 * norm(x)**2 ) otherwise. |
---|
66 | xnorm = norm(x). |
---|
67 | var estimates diag( inv(A'A) ). Omitted in this special version. |
---|
68 | outfo is a structure special to pdco.m, returning information |
---|
69 | about whether atol had to be reduced. |
---|
70 | |
---|
71 | Other potential output parameters: |
---|
72 | anorm, acond, arnorm, xnorm |
---|
73 | @endverbatim |
---|
74 | */ |
---|
75 | class ClpLsqr { |
---|
76 | private: |
---|
77 | /**@name Private member data */ |
---|
78 | //@{ |
---|
79 | //@} |
---|
80 | |
---|
81 | public: |
---|
82 | /**@name Public member data */ |
---|
83 | //@{ |
---|
84 | /// Row dimension of matrix |
---|
85 | int nrows_; |
---|
86 | /// Column dimension of matrix |
---|
87 | int ncols_; |
---|
88 | /// Pointer to Model object for this instance |
---|
89 | ClpInterior *model_; |
---|
90 | /// Diagonal array 1 |
---|
91 | double *diag1_; |
---|
92 | /// Constant diagonal 2 |
---|
93 | double diag2_; |
---|
94 | //@} |
---|
95 | |
---|
96 | /**@name Constructors and destructors */ |
---|
97 | /** Default constructor */ |
---|
98 | ClpLsqr(); |
---|
99 | |
---|
100 | /** Constructor for use with Pdco model (note modified for pdco!!!!) */ |
---|
101 | ClpLsqr(ClpInterior *model); |
---|
102 | /// Copy constructor |
---|
103 | ClpLsqr(const ClpLsqr &); |
---|
104 | /// Assignment operator. This copies the data |
---|
105 | ClpLsqr &operator=(const ClpLsqr &rhs); |
---|
106 | /** Destructor */ |
---|
107 | ~ClpLsqr(); |
---|
108 | //@} |
---|
109 | |
---|
110 | /**@name Methods */ |
---|
111 | //@{ |
---|
112 | /// Set an int parameter |
---|
113 | bool setParam(char *parmName, int parmValue); |
---|
114 | /// Call the Lsqr algorithm |
---|
115 | void do_lsqr(CoinDenseVector< double > &b, |
---|
116 | double damp, double atol, double btol, double conlim, int itnlim, |
---|
117 | bool show, Info info, CoinDenseVector< double > &x, int *istop, |
---|
118 | int *itn, Outfo *outfo, bool precon, CoinDenseVector< double > &Pr); |
---|
119 | /// Matrix-vector multiply - implemented by user |
---|
120 | void matVecMult(int, CoinDenseVector< double > *, CoinDenseVector< double > *); |
---|
121 | |
---|
122 | void matVecMult(int, CoinDenseVector< double > &, CoinDenseVector< double > &); |
---|
123 | /// diag1 - we just borrow as it is part of a CoinDenseVector<double> |
---|
124 | void borrowDiag1(double *array) |
---|
125 | { |
---|
126 | diag1_ = array; |
---|
127 | }; |
---|
128 | //@} |
---|
129 | }; |
---|
130 | #endif |
---|
131 | |
---|
132 | /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 |
---|
133 | */ |
---|