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