complex_matrices.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// This header file contains classes and inline function definitions for
27// matrices of complex numbers and their derived types
28
29// Include guards to prevent multiple inclusion of the header
30#ifndef OOMPH_COMPLEX_MATRICES_HEADER
31#define OOMPH_COMPLEX_MATRICES_HEADER
32
33// Config header generated by autoconfig
34#ifdef HAVE_CONFIG_H
35#include <oomph-lib-config.h>
36#endif
37
38#ifdef OOMPH_HAS_MPI
39#include "mpi.h"
40#endif
41
42// Of course we need the complex type
43#include <complex>
44
45// Also need the standard matrices header
46#include "matrices.h"
47
48namespace oomph
49{
50 //===================================================================
51 /// Abstract base class for matrices of complex doubles -- adds
52 /// abstract interfaces for solving, LU decomposition and
53 /// multiplication by vectors.
54 //===================================================================
56 {
57 public:
58 /// (Empty) constructor.
60
61 /// Broken copy constructor
62 ComplexMatrixBase(const ComplexMatrixBase& matrix) = delete;
63
64 /// Broken assignment operator
65 void operator=(const ComplexMatrixBase&) = delete;
66
67 /// Return the number of rows of the matrix
68 virtual unsigned long nrow() const = 0;
69
70 /// Return the number of columns of the matrix
71 virtual unsigned long ncol() const = 0;
72
73 /// virtual (empty) destructor
74 virtual ~ComplexMatrixBase() {}
75
76 /// Round brackets to give access as a(i,j) for read only
77 /// (we're not providing a general interface for component-wise write
78 /// access since not all matrix formats allow efficient direct access!)
79 virtual std::complex<double> operator()(const unsigned long& i,
80 const unsigned long& j) const = 0;
81
82 /// LU decomposition of the matrix using the appropriate
83 /// linear solver. Return the sign of the determinant
84 virtual int ludecompose()
85 {
86 throw OomphLibError(
87 "ludecompose() has not been written for this matrix class\n",
88 OOMPH_CURRENT_FUNCTION,
89 OOMPH_EXCEPTION_LOCATION);
90
91 /// Dummy return
92 return 1;
93 }
94
95 /// LU backsubstitue a previously LU-decomposed matrix;
96 /// The passed rhs will be over-written with the solution vector
97 virtual void lubksub(Vector<std::complex<double>>& rhs)
98 {
99 throw OomphLibError(
100 "lubksub() has not been written for this matrix class\n",
101 OOMPH_CURRENT_FUNCTION,
102 OOMPH_EXCEPTION_LOCATION);
103 }
104
105
106 /// Complete LU solve (replaces matrix by its LU decomposition
107 /// and overwrites RHS with solution). The default should not need
108 /// to be over-written
109 virtual void solve(Vector<std::complex<double>>& rhs);
110
111 /// Complete LU solve (Nothing gets overwritten!). The default should
112 /// not need to be overwritten
113 virtual void solve(const Vector<std::complex<double>>& rhs,
114 Vector<std::complex<double>>& soln);
115
116 /// Find the residual, i.e. r=b-Ax the residual
117 virtual void residual(const Vector<std::complex<double>>& x,
118 const Vector<std::complex<double>>& b,
119 Vector<std::complex<double>>& residual) = 0;
120
121 /// Find the maximum residual r=b-Ax -- generic version, can be
122 /// overloaded for specific derived classes where the
123 /// max. can be determined "on the fly"
124 virtual double max_residual(const Vector<std::complex<double>>& x,
125 const Vector<std::complex<double>>& rhs)
126 {
127 unsigned long n = rhs.size();
129 residual(x, rhs, res);
130 double ans = 0.0;
131 for (unsigned long i = 0; i < n; i++)
132 {
133 ans = std::max(ans, std::abs(res[i]));
134 }
135 return ans;
136 }
137
138 /// Multiply the matrix by the vector x: soln=Ax.
139 virtual void multiply(const Vector<std::complex<double>>& x,
140 Vector<std::complex<double>>& soln) = 0;
141
142 /// Multiply the transposed matrix by the vector x: soln=A^T x
143 virtual void multiply_transpose(const Vector<std::complex<double>>& x,
144 Vector<std::complex<double>>& soln) = 0;
145 };
146
147
148 //=================================================================
149 /// Class of matrices containing double complex, and stored as a
150 /// DenseMatrix<complex<double> >, but with solving functionality inherited
151 /// from the abstract ComplexMatrix class.
152 //=================================================================
154 public DenseMatrix<std::complex<double>>
155 {
156 private:
157 /// Pointer to storage for the index of permutations in the LU solve
159
160 /// Pointer to storage for the LU decomposition
161 std::complex<double>* LU_factors;
162
163 /// Boolean flag used to decided whether the LU decomposition is stored
164 /// separately, or not
166
167 /// Function that deletes the storage for the LU_factors, if it is
168 /// not the same as the matrix storage
169 void delete_lu_factors();
170
171 public:
172 /// Empty constructor, simply assign the lengths N and M to 0
174 : DenseMatrix<std::complex<double>>(),
175 Index(0),
176 LU_factors(0),
178 {
179 }
180
181 /// Constructor to build a square n by n matrix.
182 DenseComplexMatrix(const unsigned long& n)
183 : DenseMatrix<std::complex<double>>(n),
184 Index(0),
185 LU_factors(0),
187 {
188 }
189
190 /// Constructor to build a matrix with n rows and m columns.
191 DenseComplexMatrix(const unsigned long& n, const unsigned long& m)
192 : DenseMatrix<std::complex<double>>(n, m),
193 Index(0),
194 LU_factors(0),
196 {
197 }
198
199 /// Constructor to build a matrix with n rows and m columns,
200 /// with initial value initial_val
201 DenseComplexMatrix(const unsigned long& n,
202 const unsigned long& m,
203 const std::complex<double>& initial_val)
204 : DenseMatrix<std::complex<double>>(n, m, initial_val),
205 Index(0),
206 LU_factors(0),
208 {
209 }
210
211
212 /// Broken copy constructor
213 DenseComplexMatrix(const DenseComplexMatrix& matrix) = delete;
214
215 /// Broken assignment operator
216 void operator=(const DenseComplexMatrix&) = delete;
217
218 /// Return the number of rows of the matrix
219 inline unsigned long nrow() const
220 {
222 }
223
224 /// Return the number of columns of the matrix
225 inline unsigned long ncol() const
226 {
228 }
229
230 /// Overload the const version of the round-bracket access operator
231 /// for read-only access.
232 inline std::complex<double> operator()(const unsigned long& i,
233 const unsigned long& j) const
234 {
236 }
237
238 /// Overload the non-const version of the round-bracket access
239 /// operator for read-write access
240 inline std::complex<double>& operator()(const unsigned long& i,
241 const unsigned long& j)
242 {
244 }
245
246 /// Destructor, delete the storage for Index vector and LU storage (if any)
247 virtual ~DenseComplexMatrix();
248
249 /// Overload the LU decomposition.
250 /// For this storage scheme the matrix storage will be over-writeen
251 /// by its LU decomposition
252 int ludecompose();
253
254 /// Overload the LU back substitution. Note that the rhs will be
255 /// overwritten with the solution vector
256 void lubksub(Vector<std::complex<double>>& rhs);
257
258 /// Find the residual of Ax=b, ie r=b-Ax for the
259 /// "solution" x.
260 void residual(const Vector<std::complex<double>>& x,
261 const Vector<std::complex<double>>& rhs,
262 Vector<std::complex<double>>& residual);
263
264 /// Multiply the matrix by the vector x: soln=Ax
265 void multiply(const Vector<std::complex<double>>& x,
266 Vector<std::complex<double>>& soln);
267
268 /// Multiply the transposed matrix by the vector x: soln=A^T x
269 void multiply_transpose(const Vector<std::complex<double>>& x,
270 Vector<std::complex<double>>& soln);
271 };
272
273
274 //=================================================================
275 /// A class for compressed row matrices
276 //=================================================================
277 class CRComplexMatrix : public CRMatrix<std::complex<double>>,
278 public ComplexMatrixBase
279 {
280 private:
281 /// Storage for the LU factors as required by SuperLU
283
284 /// Info flag for the SuperLU solver
285 int Info;
286
287 public:
288 /// Default constructor
289 CRComplexMatrix() : CRMatrix<std::complex<double>>(), F_factors(0), Info(0)
290 {
291 // By default SuperLU solves linear systems quietly
293 }
294
295 /// Constructor: Pass vector of values, vector of column indices,
296 /// vector of row starts and number of columns (can be suppressed
297 /// for square matrices)
298 CRComplexMatrix(const Vector<std::complex<double>>& value,
300 const Vector<int>& row_start,
301 const unsigned long& n,
302 const unsigned long& m)
303 : CRMatrix<std::complex<double>>(value, column_index, row_start, n, m),
304 F_factors(0),
305 Info(0)
306 {
307 // By default SuperLU solves linear systems quietly
309 }
310
311
312 /// Broken copy constructor
313 CRComplexMatrix(const CRComplexMatrix& matrix) = delete;
314
315 /// Broken assignment operator
316 void operator=(const CRComplexMatrix&) = delete;
317
318 /// Destructor: Kill the LU decomposition if it has been computed
320 {
322 }
323
324 /// Set flag to indicate that stats are to be displayed during
325 /// solution of linear system with SuperLU
327 {
329 }
330
331 // Set flag to indicate that stats are not to be displayed during
332 /// the solve
334 {
336 }
337
338 /// Return the number of rows of the matrix
339 inline unsigned long nrow() const
340 {
342 }
343
344 /// Return the number of columns of the matrix
345 inline unsigned long ncol() const
346 {
348 }
349
350 /// Overload the round-bracket access operator for read-only access
351 inline std::complex<double> operator()(const unsigned long& i,
352 const unsigned long& j) const
353 {
355 }
356
357 /// LU decomposition using SuperLU
358 int ludecompose();
359
360 /// LU back solve for given RHS
361 void lubksub(Vector<std::complex<double>>& rhs);
362
363 /// LU clean up (perhaps this should happen in the destructor)
364 void clean_up_memory();
365
366 /// Find the residual to x of Ax=b, i.e. r=b-Ax
367 void residual(const Vector<std::complex<double>>& x,
368 const Vector<std::complex<double>>& b,
369 Vector<std::complex<double>>& residual);
370
371 /// Multiply the matrix by the vector x: soln=Ax
372 void multiply(const Vector<std::complex<double>>& x,
373 Vector<std::complex<double>>& soln);
374
375
376 /// Multiply the transposed matrix by the vector x: soln=A^T x
377 void multiply_transpose(const Vector<std::complex<double>>& x,
378 Vector<std::complex<double>>& soln);
379
380 protected:
381 /// Flag to indicate if stats are to be displayed during
382 /// solution of linear system with SuperLU
384 };
385
386
387 //=================================================================
388 /// A class for compressed column matrices that store doubles
389 //=================================================================
391 public CCMatrix<std::complex<double>>
392 {
393 private:
394 /// Storage for the LU factors as required by SuperLU
396
397 /// Info flag for the SuperLU solver
398 int Info;
399
400 protected:
401 /// Flag to indicate if stats are to be displayed during
402 /// solution of linear system with SuperLU
404
405 public:
406 /// Default constructor
407 CCComplexMatrix() : CCMatrix<std::complex<double>>(), F_factors(0), Info(0)
408 {
409 // By default SuperLU solves linear systems quietly
411 }
412
413 /// Constructor: Pass vector of values, vector of row indices,
414 /// vector of column starts and number of rows (can be suppressed
415 /// for square matrices). Number of nonzero entries is read
416 /// off from value, so make sure the vector has been shrunk
417 /// to its correct length.
418 CCComplexMatrix(const Vector<std::complex<double>>& value,
419 const Vector<int>& row_index,
421 const unsigned long& n,
422 const unsigned long& m)
423 : CCMatrix<std::complex<double>>(value, row_index, column_start, n, m),
424 F_factors(0),
425 Info(0)
426 {
427 // By default SuperLU solves linear systems quietly
429 }
430
431 /// Broken copy constructor
432 CCComplexMatrix(const CCComplexMatrix& matrix) = delete;
433
434 /// Broken assignment operator
435 void operator=(const CCComplexMatrix&) = delete;
436
437 /// Destructor: Kill the LU factors if they have been setup.
439 {
441 }
442
443 /// Set flag to indicate that stats are to be displayed during
444 /// solution of linear system with SuperLU
446 {
448 }
449
450 // Set flag to indicate that stats are not to be displayed during
451 /// the solve
453 {
455 }
456
457 /// Return the number of rows of the matrix
458 inline unsigned long nrow() const
459 {
461 }
462
463 /// Return the number of columns of the matrix
464 inline unsigned long ncol() const
465 {
467 }
468
469 /// Overload the round-bracket access operator to provide
470 /// read-only (const) access to the data
471 inline std::complex<double> operator()(const unsigned long& i,
472 const unsigned long& j) const
473 {
475 }
476
477 /// LU decomposition using SuperLU
478 int ludecompose();
479
480 /// LU back solve for given RHS
481 void lubksub(Vector<std::complex<double>>& rhs);
482
483 /// LU clean up (perhaps this should happen in the destructor)
484 void clean_up_memory();
485
486 /// Find the residulal to x of Ax=b, ie r=b-Ax
487 void residual(const Vector<std::complex<double>>& x,
488 const Vector<std::complex<double>>& b,
489 Vector<std::complex<double>>& residual);
490
491
492 /// Multiply the matrix by the vector x: soln=Ax
493 void multiply(const Vector<std::complex<double>>& x,
494 Vector<std::complex<double>>& soln);
495
496
497 /// Multiply the transposed matrix by the vector x: soln=A^T x
498 void multiply_transpose(const Vector<std::complex<double>>& x,
499 Vector<std::complex<double>>& soln);
500 };
501} // namespace oomph
502
503#endif
cstr elem_len * i
Definition: cfortran.h:603
A class for compressed column matrices that store doubles.
virtual ~CCComplexMatrix()
Destructor: Kill the LU factors if they have been setup.
int Info
Info flag for the SuperLU solver.
CCComplexMatrix(const Vector< std::complex< double > > &value, const Vector< int > &row_index, const Vector< int > &column_start, const unsigned long &n, const unsigned long &m)
Constructor: Pass vector of values, vector of row indices, vector of column starts and number of rows...
unsigned long nrow() const
Return the number of rows of the matrix.
CCComplexMatrix()
Default constructor.
unsigned long ncol() const
Return the number of columns of the matrix.
void multiply(const Vector< std::complex< double > > &x, Vector< std::complex< double > > &soln)
Multiply the matrix by the vector x: soln=Ax.
int ludecompose()
LU decomposition using SuperLU.
void multiply_transpose(const Vector< std::complex< double > > &x, Vector< std::complex< double > > &soln)
Multiply the transposed matrix by the vector x: soln=A^T x.
void clean_up_memory()
LU clean up (perhaps this should happen in the destructor)
bool Doc_stats_during_solve
Flag to indicate if stats are to be displayed during solution of linear system with SuperLU.
void operator=(const CCComplexMatrix &)=delete
Broken assignment operator.
CCComplexMatrix(const CCComplexMatrix &matrix)=delete
Broken copy constructor.
void * F_factors
Storage for the LU factors as required by SuperLU.
void disable_doc_stats()
the solve
void lubksub(Vector< std::complex< double > > &rhs)
LU back solve for given RHS.
void residual(const Vector< std::complex< double > > &x, const Vector< std::complex< double > > &b, Vector< std::complex< double > > &residual)
Find the residulal to x of Ax=b, ie r=b-Ax.
std::complex< double > operator()(const unsigned long &i, const unsigned long &j) const
Overload the round-bracket access operator to provide read-only (const) access to the data.
void enable_doc_stats()
Set flag to indicate that stats are to be displayed during solution of linear system with SuperLU.
/////////////////////////////////////////////////////////////// /////////////////////////////////////...
Definition: matrices.h:2585
std::complex< double > get_entry(const unsigned long &i, const unsigned long &j) const
Access function that will be called by the read-only round-bracket operator (const)
Definition: matrices.h:2654
int * column_start()
Access to C-style column_start array.
Definition: matrices.h:2692
int * row_index()
Access to C-style row index array.
Definition: matrices.h:2704
A class for compressed row matrices.
unsigned long nrow() const
Return the number of rows of the matrix.
void * F_factors
Storage for the LU factors as required by SuperLU.
void enable_doc_stats()
Set flag to indicate that stats are to be displayed during solution of linear system with SuperLU.
CRComplexMatrix(const Vector< std::complex< double > > &value, const Vector< int > &column_index, const Vector< int > &row_start, const unsigned long &n, const unsigned long &m)
Constructor: Pass vector of values, vector of column indices, vector of row starts and number of colu...
virtual ~CRComplexMatrix()
Destructor: Kill the LU decomposition if it has been computed.
unsigned long ncol() const
Return the number of columns of the matrix.
std::complex< double > operator()(const unsigned long &i, const unsigned long &j) const
Overload the round-bracket access operator for read-only access.
bool Doc_stats_during_solve
Flag to indicate if stats are to be displayed during solution of linear system with SuperLU.
void multiply_transpose(const Vector< std::complex< double > > &x, Vector< std::complex< double > > &soln)
Multiply the transposed matrix by the vector x: soln=A^T x.
void multiply(const Vector< std::complex< double > > &x, Vector< std::complex< double > > &soln)
Multiply the matrix by the vector x: soln=Ax.
int ludecompose()
LU decomposition using SuperLU.
CRComplexMatrix(const CRComplexMatrix &matrix)=delete
Broken copy constructor.
void disable_doc_stats()
the solve
CRComplexMatrix()
Default constructor.
void lubksub(Vector< std::complex< double > > &rhs)
LU back solve for given RHS.
void operator=(const CRComplexMatrix &)=delete
Broken assignment operator.
int Info
Info flag for the SuperLU solver.
void residual(const Vector< std::complex< double > > &x, const Vector< std::complex< double > > &b, Vector< std::complex< double > > &residual)
Find the residual to x of Ax=b, i.e. r=b-Ax.
void clean_up_memory()
LU clean up (perhaps this should happen in the destructor)
A class for compressed row matrices, a sparse storage format Once again the recursive template trick ...
Definition: matrices.h:682
int * column_index()
Access to C-style column index array.
Definition: matrices.h:797
int * row_start()
Access to C-style row_start array.
Definition: matrices.h:785
std::complex< double > get_entry(const unsigned long &i, const unsigned long &j) const
Access function that will be called by the read-only round-bracket operator (const)
Definition: matrices.h:748
Abstract base class for matrices of complex doubles – adds abstract interfaces for solving,...
virtual int ludecompose()
LU decomposition of the matrix using the appropriate linear solver. Return the sign of the determinan...
virtual void solve(Vector< std::complex< double > > &rhs)
Complete LU solve (replaces matrix by its LU decomposition and overwrites RHS with solution)....
virtual std::complex< double > operator()(const unsigned long &i, const unsigned long &j) const =0
Round brackets to give access as a(i,j) for read only (we're not providing a general interface for co...
virtual ~ComplexMatrixBase()
virtual (empty) destructor
virtual void multiply(const Vector< std::complex< double > > &x, Vector< std::complex< double > > &soln)=0
Multiply the matrix by the vector x: soln=Ax.
virtual unsigned long ncol() const =0
Return the number of columns of the matrix.
virtual void multiply_transpose(const Vector< std::complex< double > > &x, Vector< std::complex< double > > &soln)=0
Multiply the transposed matrix by the vector x: soln=A^T x.
void operator=(const ComplexMatrixBase &)=delete
Broken assignment operator.
virtual void lubksub(Vector< std::complex< double > > &rhs)
LU backsubstitue a previously LU-decomposed matrix; The passed rhs will be over-written with the solu...
virtual double max_residual(const Vector< std::complex< double > > &x, const Vector< std::complex< double > > &rhs)
Find the maximum residual r=b-Ax – generic version, can be overloaded for specific derived classes wh...
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
virtual void residual(const Vector< std::complex< double > > &x, const Vector< std::complex< double > > &b, Vector< std::complex< double > > &residual)=0
Find the residual, i.e. r=b-Ax the residual.
ComplexMatrixBase()
(Empty) constructor.
ComplexMatrixBase(const ComplexMatrixBase &matrix)=delete
Broken copy constructor.
Class of matrices containing double complex, and stored as a DenseMatrix<complex<double> >,...
bool Overwrite_matrix_storage
Boolean flag used to decided whether the LU decomposition is stored separately, or not.
virtual ~DenseComplexMatrix()
Destructor, delete the storage for Index vector and LU storage (if any)
unsigned long ncol() const
Return the number of columns of the matrix.
Vector< long > * Index
Pointer to storage for the index of permutations in the LU solve.
void residual(const Vector< std::complex< double > > &x, const Vector< std::complex< double > > &rhs, Vector< std::complex< double > > &residual)
Find the residual of Ax=b, ie r=b-Ax for the "solution" x.
std::complex< double > & operator()(const unsigned long &i, const unsigned long &j)
Overload the non-const version of the round-bracket access operator for read-write access.
int ludecompose()
Overload the LU decomposition. For this storage scheme the matrix storage will be over-writeen by its...
DenseComplexMatrix(const unsigned long &n, const unsigned long &m, const std::complex< double > &initial_val)
Constructor to build a matrix with n rows and m columns, with initial value initial_val.
void multiply_transpose(const Vector< std::complex< double > > &x, Vector< std::complex< double > > &soln)
Multiply the transposed matrix by the vector x: soln=A^T x.
DenseComplexMatrix(const DenseComplexMatrix &matrix)=delete
Broken copy constructor.
void delete_lu_factors()
Function that deletes the storage for the LU_factors, if it is not the same as the matrix storage.
DenseComplexMatrix(const unsigned long &n, const unsigned long &m)
Constructor to build a matrix with n rows and m columns.
void lubksub(Vector< std::complex< double > > &rhs)
Overload the LU back substitution. Note that the rhs will be overwritten with the solution vector.
void operator=(const DenseComplexMatrix &)=delete
Broken assignment operator.
DenseComplexMatrix(const unsigned long &n)
Constructor to build a square n by n matrix.
std::complex< double > operator()(const unsigned long &i, const unsigned long &j) const
Overload the const version of the round-bracket access operator for read-only access.
DenseComplexMatrix()
Empty constructor, simply assign the lengths N and M to 0.
std::complex< double > * LU_factors
Pointer to storage for the LU decomposition.
unsigned long nrow() const
Return the number of rows of the matrix.
void multiply(const Vector< std::complex< double > > &x, Vector< std::complex< double > > &soln)
Multiply the matrix by the vector x: soln=Ax.
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
Definition: matrices.h:386
std::complex< double > get_entry(const unsigned long &i, const unsigned long &j) const
The access function the will be called by the read-only (const version) round-bracket operator.
Definition: matrices.h:457
std::complex< double > & entry(const unsigned long &i, const unsigned long &j)
The access function that will be called by the read-write round-bracket operator.
Definition: matrices.h:447
An OomphLibError object which should be thrown when an run-time error is encountered....
T * value()
Access to C-style value array.
Definition: matrices.h:616
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...