source: trunk/Clp/src/ClpLsqr.hpp @ 1304

Last change on this file since 1304 was 1171, checked in by forrest, 12 years ago

put back pdco in case wanted

File size: 4.4 KB
Line 
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/**
12This 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 */
73class ClpLsqr{
74private:
75   /**@name Private member data */
76   //@{
77   //@}
78
79public:
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   
Note: See TracBrowser for help on using the repository browser.