Chombo + EB  3.2
LAPACKMatrix.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
8  */
9 #endif
10 
11 #ifndef _LAPACKMATRIX_H_
12 #define _LAPACKMATRIX_H_
13 
14 #include "REAL.H"
15 #include "LoHiSide.H"
16 #include <utility>
17 #include "BaseNamespaceHeader.H"
18 
19 ///
20 /**
21  Getting sick of writing the basics here over and over.
22  Silly class but it will cut down on the typing.
23  */
25 {
26 public:
27 
28 
29  /// main constructor. Matrix of Chombo Real type. values are unintitialized
30  LAPACKMatrix(int a_nrow, int a_ncol):m_nrow(a_nrow),m_ncol(a_ncol),m_data(new Real[a_nrow*a_ncol]){}
31 
32  ///
33  void define(int nrow, int ncol)
34  {
35  m_nrow = nrow;
36  m_ncol = ncol;
37  if(m_data != nullptr)
38  {
39  delete m_data;
40  }
41  m_data =new Real[m_nrow*m_ncol];
42  }
43 
44 
45  /// alias constructor. Use a_data as m_data, does not delete it in destructor
46  LAPACKMatrix( int a_nrow, int a_ncol, Real* a_data)
47  :m_nrow(a_nrow),m_ncol(a_ncol),m_data(a_data),m_alias(true) {;}
48 
49  /// null constructor builds matrix with nullptr and size=[0,0]
50  LAPACKMatrix():m_nrow(0),m_ncol(0),m_data(nullptr){}
51 
52  /// deep copy constructor
53  LAPACKMatrix(const LAPACKMatrix& a_input);
54 
55  /// move copy constructor
57  :m_nrow(a_input.m_nrow),m_ncol(a_input.m_ncol),m_data(a_input.m_data),m_alias(a_input.m_alias)
58  {a_input.m_data=nullptr;}
59 
60 
61  void clear()
62  {
63  define(0, 0);
64  }
65 
66  ///return sqrt(sum of squares of all values)
67  Real normLTwo() const;
68 
69  ///
70  /**
71  if(i==j) M=1 else M=0
72  */
73  void setToIdentity();
74 
75  ///
76  ~LAPACKMatrix();
77 
78  ///
79  void setVal(const Real& a_val);
80 
81  Real* dataPtr();
82 
83  ///
84  const Real* dataPtr() const;
85 
86  ///
87  /**
88  make nrows = a_nrows etc. Useful when you want to discard the last few rows of the matrix.
89  */
90  void truncate(int a_nrow,int a_ncol);
91 
92 
93  ///
94  /**
95  Get the maximum absolute value of the matrix (for iterative solves)
96  */
97  Real maxNorm() const;
98 
99  ///
100  const Real& operator() (int irow, int icol) const;
101 
102  ///
103  Real & operator() (int irow, int icol);
104 
105  //return row, col
106  std::pair<int, int> dims() const
107  {
108  std::pair<int, int> retval;
109  retval.first = m_nrow;
110  retval.second= m_ncol;
111  return retval;
112  }
113 
114  /// deep assign.
115  LAPACKMatrix& operator=(const LAPACKMatrix& a_matrix);
116 
117  /// move assignment
118  inline LAPACKMatrix& operator=(LAPACKMatrix&& a_matrix)
119  {
120  if(&a_matrix != this)
121  {
122  m_nrow=a_matrix.m_nrow;
123  m_ncol=a_matrix.m_ncol;
124  std::swap(m_data,a_matrix.m_data);
125  std::swap(m_alias,a_matrix.m_alias);
126  }
127  return *this;
128  }
129 
130 
131  ///
132  LAPACKMatrix& operator+=(const LAPACKMatrix& a_matrix);
133 
134  ///
135  LAPACKMatrix& operator-=(const LAPACKMatrix& a_matrix);
136 
137  ///
138  LAPACKMatrix& operator*=(const Real& a_scalingFactor);
139 
140 
141  ///
142  int offset(int irow, int icol) const;
143 
144  ///
145  void poutAll() const;
146 
147  ///
148  void poutDiag() const;
149 
150  ///
151  void poutDiagMatlab() const;
152 
153  ///
154  void poutMatlab() const;
155 
156  ///inverts this matix
157  /**
158  fails if matrix is not square
159  if return value != 0, probably a singular matrix
160  */
161  int invert();
162 
163  /// inverts this matrix using SVD
164  /**
165  * Get Ainverse using least squares with svd
166  * if return value != 0, lapack was unhappy somehow
167  */
168  int invertUsingSVD(int a_maxiter, Real a_tol);
169 
170  ///
171  //pseudoinverse = (AT A)^(-1) AT
172  int pseudoInvertUsingSVD(int a_maxiter, Real a_tol);
173 
174  // pseudoinverse: X such that R*X = Q^T, where A = Q*R; hence A*X = Q*R*X = Q*Q^T = I
175  int pseudoInvertUsingQR();
176 
177  /// inverts this matrix using least squares
178  /**
179  * Get Ainverse using least squares with svd
180  * if return value != 0, lapack was unhappy somehow
181  */
183 
184 
185  //for tiny cells, set moments to their limits (1 0 0 0...)
186  void setSmallCellRow(const int& irow)
187  {
188  LAPACKMatrix& Mvol = *this;
189  Mvol(irow, 0) = 1.0;
190  for(int icol = 1; icol < m_ncol; icol++)
191  {
192  Mvol(irow, icol) = 0.0;
193  }
194  }
195 
196 
197  ///
198  void transpose();
199 
200 
201  void checkConditionNumber() const;
202 
204 
205  ///
206  /**
207  sets product = a_left* a_right
208  fails if a_left.m_col != a_right.m_rows
209  */
210  friend void multiply(LAPACKMatrix& a_product,
211  const LAPACKMatrix& a_left,
212  const LAPACKMatrix& a_right);
213 
214  ///below stuff is shamelessly stolen from lapackwrapper class
215 
216  ///
217  /**
218  Solves A*X = B using general least squares, for each column of B
219  */
221 
222  ///
223  /**
224  Solves A'*X = B using least squares, for vector b. Answer goes
225  back into B I think
226  */
228 
229  ///
230  /**
231  * Solves A^T X = B using least squares with SVD, for vector b
232  */
233  friend int solveLSTSVD(LAPACKMatrix& A, LAPACKMatrix& B, int a_maxiter, Real a_tol);
234 
235  ///
236  /**
237  * Solves A*X = B using least squares with SVD, for X
238  */
239  friend int solveLSTSVD(LAPACKMatrix& X, const LAPACKMatrix& A, const LAPACKMatrix& B, int a_maxiter, Real a_tol);
240 
241  ///
242  friend int solveLSTSVDOnce(LAPACKMatrix& X, LAPACKMatrix& B);
243 
244  ///
245  /**
246  Solves A*X = B using least squares with SVD, for X
247  */
248  friend int solveLSTSVDOnce(LAPACKMatrix& X, const LAPACKMatrix& A, const LAPACKMatrix& B);
249 
250  ///
251  /**
252  * Solves equality constrained least squares problem
253  * Find x, s.t. min norm(A x - c) with B x = d
254  */
256 
257  ///
258  /**
259  * Solves A'*X = B using reduced rank least squares, for vector b
260  */
262 
263  ///
264  /**
265  Following Lapack, gets inverse of condition number. Returning a number
266  near zero means the matrix is not really solvable.
267  */
269 
271 
272  /// turn on if you want every solve to check the condition number
274 
275  ///
276  static bool s_verbose;
277  ///
278  static bool s_outputStenData;
279 
280 private:
281 
282  int m_nrow;
283  int m_ncol;
285  bool m_alias=false;
286 };
287 ///
289 
290 ///
291 /**
292  sets product = a_left* a_right
293  fails if a_left.m_col != a_right.m_rows
294 */
295 void multiply(LAPACKMatrix& a_product,
296  const LAPACKMatrix& a_left,
297  const LAPACKMatrix& a_right);
298 
299 ///below stuff is shamelessly stolen from lapackwrapper class
300 
301 ///
302 /**
303  Solves A*X = B using general least squares, for each column of B
304 */
306 
307 ///
308 /**
309  Solves A'*X = B using least squares, for vector b
310 */
312 
313 ///
314 int solveLSTSVDOnce(LAPACKMatrix& X, const LAPACKMatrix& A, const LAPACKMatrix& B);
315 
316 ///
317 /**
318  * Solves A^T X = B using least squares with SVD, for vector b
319  */
320 int solveLSTSVD(LAPACKMatrix& A, LAPACKMatrix& B, int a_maxiter, Real a_tol);
321 
322 ///
323 int solveLSTSVDOnce(LAPACKMatrix& X, const LAPACKMatrix& A, const LAPACKMatrix& B);
324 
325 ///
326 /**
327  * Solves equality constrained least squares problem
328  * Find x, s.t. min norm(A x - c) with B x = d
329  */
331 
332 ///
333 /**
334  * Solves A'*X = B using reduced rank least squares, for vector b
335  */
337 
338 
339 #include "BaseNamespaceFooter.H"
340 #endif
void poutDiag() const
LAPACKMatrix(int a_nrow, int a_ncol)
main constructor. Matrix of Chombo Real type. values are unintitialized
Definition: LAPACKMatrix.H:30
int solveReducedRankLS(LAPACKMatrix &A, LAPACKMatrix &b)
LAPACKMatrix(LAPACKMatrix &&a_input)
move copy constructor
Definition: LAPACKMatrix.H:56
void transpose()
int solveEqualityConstrainedLS(LAPACKMatrix &A, LAPACKMatrix &c, LAPACKMatrix &B, LAPACKMatrix &d, LAPACKMatrix &x)
int invertUsingLeastSquares()
inverts this matrix using least squares
int solveLeastSquaresTranspose(LAPACKMatrix &A, LAPACKMatrix &B)
friend int solveReducedRankLS(LAPACKMatrix &A, LAPACKMatrix &b)
LAPACKMatrix & operator*=(const Real &a_scalingFactor)
friend Real getInverseOfUpperTriangularConditionNumber(const LAPACKMatrix &A)
int pseudoInvertUsingSVD(int a_maxiter, Real a_tol)
void define(int nrow, int ncol)
Definition: LAPACKMatrix.H:33
LAPACKMatrix & operator=(LAPACKMatrix &&a_matrix)
move assignment
Definition: LAPACKMatrix.H:118
int pseudoInvertUsingQR()
friend int solveLSTSVDOnce(LAPACKMatrix &X, LAPACKMatrix &B)
const Real & operator()(int irow, int icol) const
static bool s_verbose
Definition: LAPACKMatrix.H:276
bool m_alias
Definition: LAPACKMatrix.H:285
void multiply(LAPACKMatrix &a_product, const LAPACKMatrix &a_left, const LAPACKMatrix &a_right)
void const char const int const int const int const Real const Real * A
Definition: Lapack.H:83
friend int solveLeastSquares(LAPACKMatrix &A, LAPACKMatrix &B)
below stuff is shamelessly stolen from lapackwrapper class
void setVal(const Real &a_val)
friend int solveEqualityConstrainedLS(LAPACKMatrix &A, LAPACKMatrix &c, LAPACKMatrix &B, LAPACKMatrix &d, LAPACKMatrix &x)
static bool s_outputStenData
Definition: LAPACKMatrix.H:278
int offset(int irow, int icol) const
int solveLSTSVD(LAPACKMatrix &A, LAPACKMatrix &B, int a_maxiter, Real a_tol)
double Real
Definition: REAL.H:33
Real * dataPtr()
friend int solveLSTSVD(LAPACKMatrix &A, LAPACKMatrix &B, int a_maxiter, Real a_tol)
Definition: LAPACKMatrix.H:24
int m_nrow
Definition: LAPACKMatrix.H:282
Real getInverseOfConditionNumber(const LAPACKMatrix &A)
Real maxNorm() const
Real * m_data
Definition: LAPACKMatrix.H:284
LAPACKMatrix(int a_nrow, int a_ncol, Real *a_data)
alias constructor. Use a_data as m_data, does not delete it in destructor
Definition: LAPACKMatrix.H:46
int invertUsingSVD(int a_maxiter, Real a_tol)
inverts this matrix using SVD
LAPACKMatrix & operator=(const LAPACKMatrix &a_matrix)
deep assign.
LAPACKMatrix & operator+=(const LAPACKMatrix &a_matrix)
void truncate(int a_nrow, int a_ncol)
static bool s_checkConditionNumber
turn on if you want every solve to check the condition number
Definition: LAPACKMatrix.H:273
void checkConditionNumber() const
Real normLTwo() const
return sqrt(sum of squares of all values)
void poutDiagMatlab() const
void clear()
Definition: LAPACKMatrix.H:61
void setToIdentity()
void poutAll() const
int solveLSTSVDOnce(LAPACKMatrix &X, const LAPACKMatrix &A, const LAPACKMatrix &B)
int m_ncol
Definition: LAPACKMatrix.H:283
void const char const int const int const int const Real const Real const int const Real * B
Definition: Lapack.H:83
friend Real getInverseOfConditionNumber(const LAPACKMatrix &A)
std::pair< int, int > dims() const
Definition: LAPACKMatrix.H:106
void setSmallCellRow(const int &irow)
Definition: LAPACKMatrix.H:186
void checkUpperTriangularConditionNumber() const
LAPACKMatrix & operator-=(const LAPACKMatrix &a_matrix)
friend int solveLeastSquaresTranspose(LAPACKMatrix &A, LAPACKMatrix &B)
void poutMatlab() const
LAPACKMatrix()
null constructor builds matrix with nullptr and size=[0,0]
Definition: LAPACKMatrix.H:50
int solveLeastSquares(LAPACKMatrix &A, LAPACKMatrix &B)
below stuff is shamelessly stolen from lapackwrapper class
int invert()
inverts this matix
friend void multiply(LAPACKMatrix &a_product, const LAPACKMatrix &a_left, const LAPACKMatrix &a_right)