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