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-2023 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 and their derived types
28 
29 // Include guards to prevent multiple inclusion of the header
30 #ifndef OOMPH_MATRICES_HEADER
31 #define OOMPH_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 
43 // Needed for g++ in some cases
44 #include <iomanip>
45 
46 // oomph-lib headers
47 #include "Vector.h"
48 #include "oomph_utilities.h"
50 #include "double_vector.h"
51 
52 
53 #ifdef OOMPH_HAS_TRILINOS
54 #include "trilinos_helpers.h"
55 #endif
56 
57 namespace oomph
58 {
59 // Initialise dense pointer-based matrices/tensors?
60 #define OOMPH_INITIALISE_DENSE_MATRICES
61 #undef OOMPH_INITIALISE_DENSE_MATRICES
62 
63  //=================================================================
64  /// Abstract base class for matrices, templated by
65  /// the type of object that is stored in them and the type of matrix.
66  /// The MATRIX_TYPE template argument is used as part of the
67  /// Curiously Recurring Template Pattern, see
68  /// http://en.wikipedia.org/wiki/Curiously_Recurring_Template_Pattern
69  /// The pattern is used to force the inlining of the round bracket access
70  /// functions by ensuring that they are NOT virtual functions.
71  //=================================================================
72  template<class T, class MATRIX_TYPE>
73  class Matrix
74  {
75  protected:
76  /// Range check to catch when an index is out of bounds, if so, it
77  /// issues a warning message and dies by throwing an \c OomphLibError
78  void range_check(const unsigned long& i, const unsigned long& j) const
79  {
80  if (i >= nrow())
81  {
82  std::ostringstream error_message;
83  error_message << "Range Error: i=" << i << " is not in the range (0,"
84  << nrow() - 1 << ")." << std::endl;
85 
86  throw OomphLibError(error_message.str(),
87  OOMPH_CURRENT_FUNCTION,
88  OOMPH_EXCEPTION_LOCATION);
89  }
90  else if (j >= ncol())
91  {
92  std::ostringstream error_message;
93  error_message << "Range Error: j=" << j << " is not in the range (0,"
94  << ncol() - 1 << ")." << std::endl;
95 
96  throw OomphLibError(error_message.str(),
97  OOMPH_CURRENT_FUNCTION,
98  OOMPH_EXCEPTION_LOCATION);
99  }
100  }
101 
102 
103  public:
104  /// (Empty) constructor
105  Matrix() {}
106 
107  /// Broken copy constructor
108  Matrix(const Matrix& matrix) = delete;
109 
110  /// Broken assignment operator
111  void operator=(const Matrix&) = delete;
112 
113  /// Virtual (empty) destructor
114  virtual ~Matrix() {}
115 
116  /// Return the number of rows of the matrix
117  virtual unsigned long nrow() const = 0;
118 
119  /// Return the number of columns of the matrix
120  virtual unsigned long ncol() const = 0;
121 
122  /// Round brackets to give access as a(i,j) for read only
123  /// (we're not providing a general interface for component-wise write
124  /// access since not all matrix formats allow efficient direct access!)
125  /// The function uses the MATRIX_TYPE template parameter to call the
126  /// get_entry() function which must be defined in all derived classes
127  /// that are to be fully instantiated.
128  inline T operator()(const unsigned long& i, const unsigned long& j) const
129  {
130  return static_cast<MATRIX_TYPE const*>(this)->get_entry(i, j);
131  }
132 
133  /// Round brackets to give access as a(i,j) for read-write
134  /// access.
135  /// The function uses the MATRIX_TYPE template parameter to call the
136  /// entry() function which must be defined in all derived classes
137  /// that are to be fully instantiated. If the particular Matrix does
138  /// not allow write access, the function should break with an error
139  /// message.
140  inline T& operator()(const unsigned long& i, const unsigned long& j)
141  {
142  return static_cast<MATRIX_TYPE*>(this)->entry(i, j);
143  }
144 
145  /// Output function to print a matrix row-by-row, in the form
146  /// a(0,0) a(0,1) ...
147  /// a(1,0) a(1,1) ...
148  /// ...
149  /// to the stream outfile.
150  /// Broken virtual since it might not be sensible to implement this for
151  /// some sparse matrices.
152  virtual void output(std::ostream& outfile) const
153  {
154  throw OomphLibError(
155  "Output function is not implemented for this matrix class",
156  OOMPH_CURRENT_FUNCTION,
157  OOMPH_EXCEPTION_LOCATION);
158  }
159 
160  /// Output the "bottom right" entry regardless of it being
161  /// zero or not (this allows automatic detection of matrix size in
162  /// e.g. matlab, python).
163  /// This functionality was moved from the function
164  /// sparse_indexed_output(...) because at the moment, generalisation of
165  /// this functionality does not work in parallel. CRDoubleMatrix has an
166  /// nrow() function but it should it should use nrow_local() - which is the
167  /// N variable in the underlaying CRMatrix.
169  std::ostream& outfile) const = 0;
170 
171  /// Indexed output function to print a matrix to the stream outfile
172  /// as i,j,a(i,j) for a(i,j)!=0 only.
173  virtual void sparse_indexed_output_helper(std::ostream& outfile) const = 0;
174 
175 
176  /// Indexed output function to print a matrix to the stream outfile
177  /// as i,j,a(i,j) for a(i,j)!=0 only with specified precision (if
178  /// precision=0 then nothing is changed). If optional boolean flag is set
179  /// to true we also output the "bottom right" entry regardless of it being
180  /// zero or not (this allows automatic detection of matrix size in
181  /// e.g. matlab, python).
183  std::ostream& outfile,
184  const unsigned& precision = 0,
185  const bool& output_bottom_right_zero = false) const
186  {
187  // Implemented as a wrapper around "sparse_indexed_output(std::ostream)"
188  // so that only one output helper function is needed in derived classes.
189 
190  // We can't have separate functions for only "output_bottom_right_zero"
191  // because people often write false as "0" and then C++ would pick the
192  // wrong function.
193 
194  // If requested set the new precision and store the previous value.
195  unsigned old_precision = 0;
196  if (precision != 0)
197  {
198  old_precision = outfile.precision();
199  outfile.precision(precision);
200  }
201 
202  // Output as normal using the helper function defined in each matrix class
204 
205  // If requested and there is no output for the last entry then output a
206  // zero entry.
207  if (output_bottom_right_zero && ncol() > 0 && nrow() > 0)
208  {
209  // Output as normal using the helper function defined
210  // in each matrix class
212  }
213 
214  // Restore the old value of the precision if we changed it
215  if (precision != 0)
216  {
217  outfile.precision(old_precision);
218  }
219  }
220 
221  /// Indexed output function to print a matrix to the file named
222  /// filename as i,j,a(i,j) for a(i,j)!=0 only with specified precision. If
223  /// optional boolean flag is set to true we also output the "bottom right"
224  /// entry regardless of it being zero or not (this allows automatic
225  /// detection of matrix size in e.g. matlab, python).
227  std::string filename,
228  const unsigned& precision = 0,
229  const bool& output_bottom_right_zero = false) const
230  {
231  // Implemented as a wrapper around "sparse_indexed_output(std::ostream)"
232  // so that only one output function needs to be written in matrix
233  // subclasses.
234 
235  // Open file
236  std::ofstream some_file(filename.c_str());
237 
238  // Output as normal
239  sparse_indexed_output(some_file, precision, output_bottom_right_zero);
240 
241  // Close file
242  some_file.close();
243  }
244  };
245 
246 
247  /// //////////////////////////////////////////////////////////////////
248  /// //////////////////////////////////////////////////////////////////
249  /// //////////////////////////////////////////////////////////////////
250 
251 
252  // Forward definition of the linear solver class
253  class LinearSolver;
254 
255  //=============================================================================
256  /// Abstract base class for matrices of doubles -- adds
257  /// abstract interfaces for solving, LU decomposition and
258  /// multiplication by vectors.
259  //=============================================================================
261  {
262  protected:
263  // Pointer to a linear solver
265 
266  // Pointer to a default linear solver
268 
269  public:
270  /// (Empty) constructor.
272 
273  /// Broken copy constructor
274  DoubleMatrixBase(const DoubleMatrixBase& matrix) = delete;
275 
276  /// Broken assignment operator
277  void operator=(const DoubleMatrixBase&) = delete;
278 
279  /// Return the number of rows of the matrix
280  virtual unsigned long nrow() const = 0;
281 
282  /// Return the number of columns of the matrix
283  virtual unsigned long ncol() const = 0;
284 
285  /// virtual (empty) destructor
286  virtual ~DoubleMatrixBase() {}
287 
288  /// Round brackets to give access as a(i,j) for read only
289  /// (we're not providing a general interface for component-wise write
290  /// access since not all matrix formats allow efficient direct access!)
291  virtual double operator()(const unsigned long& i,
292  const unsigned long& j) const = 0;
293 
294 
295  /// Return a pointer to the linear solver object
297  {
298  return Linear_solver_pt;
299  }
300 
301  /// Return a pointer to the linear solver object (const version)
303  {
304  return Linear_solver_pt;
305  }
306 
307  /// Complete LU solve (replaces matrix by its LU decomposition
308  /// and overwrites RHS with solution). The default should not need
309  /// to be over-written
310  void solve(DoubleVector& rhs);
311 
312  /// Complete LU solve (Nothing gets overwritten!). The default should
313  /// not need to be overwritten
314  void solve(const DoubleVector& rhs, DoubleVector& soln);
315 
316  /// Complete LU solve (replaces matrix by its LU decomposition
317  /// and overwrites RHS with solution). The default should not need
318  /// to be over-written
319  void solve(Vector<double>& rhs);
320 
321  /// Complete LU solve (Nothing gets overwritten!). The default should
322  /// not need to be overwritten
323  void solve(const Vector<double>& rhs, Vector<double>& soln);
324 
325  /// Find the residual, i.e. r=b-Ax the residual
326  virtual void residual(const DoubleVector& x,
327  const DoubleVector& b,
328  DoubleVector& residual_)
329  {
330  // compute residual = Ax
331  this->multiply(x, residual_);
332 
333  // set residual to -residual (-Ax)
334  unsigned nrow_local = residual_.nrow_local();
335  double* residual_pt = residual_.values_pt();
336  for (unsigned i = 0; i < nrow_local; i++)
337  {
338  residual_pt[i] = -residual_pt[i];
339  }
340 
341  // set residual = b + residuals
342  residual_ += b;
343  }
344 
345  /// Find the maximum residual r=b-Ax -- generic version, can be
346  /// overloaded for specific derived classes where the
347  /// max. can be determined "on the fly"
348  virtual double max_residual(const DoubleVector& x, const DoubleVector& rhs)
349  {
350  DoubleVector res;
351  residual(x, rhs, res);
352  return res.max();
353  }
354 
355  /// Multiply the matrix by the vector x: soln=Ax.
356  virtual void multiply(const DoubleVector& x, DoubleVector& soln) const = 0;
357 
358  /// Multiply the transposed matrix by the vector x: soln=A^T x
359  virtual void multiply_transpose(const DoubleVector& x,
360  DoubleVector& soln) const = 0;
361 
362  /// For every row, find the maximum absolute value of the
363  /// entries in this row. Set all values that are less than alpha times
364  /// this maximum to zero and return the resulting matrix in
365  /// reduced_matrix. Note: Diagonal entries are retained regardless
366  /// of their size.
367  // virtual void matrix_reduction(const double &alpha,
368  // DoubleMatrixBase& reduced_matrix)=0;
369  };
370 
371 
372  /// ////////////////////////////////////////////////////////////////////////////
373  /// ////////////////////////////////////////////////////////////////////////////
374  /// ////////////////////////////////////////////////////////////////////////////
375 
376 
377  //======================================================================
378  /// Class for dense matrices, storing all the values of the
379  /// matrix as a pointer to a pointer with assorted output functions
380  /// inherited from Matrix<T>. The curious recursive template pattern is
381  /// used here to pass the specific class to the base class so that
382  /// round bracket access can be inlined.
383  //======================================================================
384  template<class T>
385  class DenseMatrix : public Matrix<T, DenseMatrix<T>>
386  {
387  protected:
388  /// Internal representation of matrix as a pointer to data
390 
391  /// Number of rows
392  unsigned long N;
393 
394  /// Number of columns
395  unsigned long M;
396 
397  public:
398  /// Empty constructor, simply assign the lengths N and M to 0
399  DenseMatrix() : Matrixdata(0), N(0), M(0) {}
400 
401  /// Copy constructor: Deep copy!
402  DenseMatrix(const DenseMatrix& source_matrix)
403  {
404  // Set row and column lengths
405  N = source_matrix.nrow();
406  M = source_matrix.ncol();
407  // Assign space for the data
408  Matrixdata = new T[N * M];
409  // Copy the data across from the other matrix
410  for (unsigned long i = 0; i < N; i++)
411  {
412  for (unsigned long j = 0; j < M; j++)
413  {
414  Matrixdata[M * i + j] = source_matrix(i, j);
415  }
416  }
417  }
418 
419  /// Copy assignment
420  DenseMatrix& operator=(const DenseMatrix& source_matrix)
421  {
422  // Don't create a new matrix if the assignment is the identity
423  if (this != &source_matrix)
424  {
425  // Check row and column length
426  unsigned long n = source_matrix.nrow();
427  unsigned long m = source_matrix.ncol();
428  if ((N != n) || (M != m))
429  {
430  resize(n, m);
431  }
432  // Copy entries across from the other matrix
433  for (unsigned long i = 0; i < N; i++)
434  {
435  for (unsigned long j = 0; j < M; j++)
436  {
437  (*this)(i, j) = source_matrix(i, j);
438  }
439  }
440  }
441  // Return reference to object itself (i.e. de-reference this pointer)
442  return *this;
443  }
444 
445  /// The access function that will be called by the read-write
446  /// round-bracket operator.
447  inline T& entry(const unsigned long& i, const unsigned long& j)
448  {
449 #ifdef RANGE_CHECKING
450  this->range_check(i, j);
451 #endif
452  return Matrixdata[M * i + j];
453  }
454 
455  /// The access function the will be called by the read-only
456  /// (const version) round-bracket operator.
457  inline T get_entry(const unsigned long& i, const unsigned long& j) const
458  {
459 #ifdef RANGE_CHECKING
460  this->range_check(i, j);
461 #endif
462  return Matrixdata[M * i + j];
463  }
464 
465  /// Constructor to build a square n by n matrix
466  DenseMatrix(const unsigned long& n);
467 
468  /// Constructor to build a matrix with n rows and m columns
469  DenseMatrix(const unsigned long& n, const unsigned long& m);
470 
471  /// Constructor to build a matrix with n rows and m columns,
472  /// with initial value initial_val
473  DenseMatrix(const unsigned long& n,
474  const unsigned long& m,
475  const T& initial_val);
476 
477  /// Destructor, clean up the matrix data
478  virtual ~DenseMatrix()
479  {
480  delete[] Matrixdata;
481  Matrixdata = 0;
482  }
483 
484  /// Return the number of rows of the matrix
485  inline unsigned long nrow() const
486  {
487  return N;
488  }
489 
490  /// Return the number of columns of the matrix
491  inline unsigned long ncol() const
492  {
493  return M;
494  }
495 
496  /// Resize to a square nxn matrix;
497  /// any values already present will be transfered
498  void resize(const unsigned long& n)
499  {
500  resize(n, n);
501  }
502 
503  /// Resize to a non-square n x m matrix;
504  /// any values already present will be transfered
505  void resize(const unsigned long& n, const unsigned long& m);
506 
507  /// Resize to a non-square n x m matrix and initialize the
508  /// new values to initial_value.
509  void resize(const unsigned long& n,
510  const unsigned long& m,
511  const T& initial_value);
512 
513  /// Initialize all values in the matrix to val.
514  void initialise(const T& val)
515  {
516  for (unsigned long i = 0; i < (N * M); ++i)
517  {
518  Matrixdata[i] = val;
519  }
520  }
521 
522  /// Output function to print a matrix row-by-row to the stream outfile
523  void output(std::ostream& outfile) const;
524 
525  /// Output function to print a matrix row-by-row to a file. Specify
526  /// filename.
527  void output(std::string filename) const;
528 
529  /// Indexed output function to print a matrix to the
530  /// stream outfile as i,j,a(i,j)
531  void indexed_output(std::ostream& outfile) const;
532 
533  /// Indexed output function to print a matrix to a
534  /// file as i,j,a(i,j). Specify filename.
535  void indexed_output(std::string filename) const;
536 
537  /// Output the "bottom right" entry regardless of it being
538  /// zero or not (this allows automatic detection of matrix size in
539  /// e.g. matlab, python).
540  void output_bottom_right_zero_helper(std::ostream& outfile) const;
541 
542  /// Indexed output function to print a matrix to the
543  /// stream outfile as i,j,a(i,j) for a(i,j)!=0 only.
544  void sparse_indexed_output_helper(std::ostream& outfile) const;
545  };
546 
547 
548  /// ////////////////////////////////////////////////////////////////
549  /// ////////////////////////////////////////////////////////////////
550  /// ////////////////////////////////////////////////////////////////
551 
552 
553  //================================================================
554  /// Class for sparse matrices, that store only the non-zero values
555  /// in a linear array in memory. The details of the array indexing
556  /// vary depending on the storage scheme used. The MATRIX_TYPE
557  /// template parameter for use in the curious recursive template
558  /// pattern is included and passed directly to the base Matrix class.
559  //=================================================================
560  template<class T, class MATRIX_TYPE>
561  class SparseMatrix : public Matrix<T, MATRIX_TYPE>
562  {
563  protected:
564  /// Internal representation of the matrix values, a pointer
566 
567  /// Number of rows
568  unsigned long N;
569 
570  /// Number of columns
571  unsigned long M;
572 
573  /// Number of non-zero values (i.e. size of Value array)
574  unsigned long Nnz;
575 
576  /// Dummy zero
577  static T Zero;
578 
579  public:
580  /// Default constructor
581  SparseMatrix() : Value(0), N(0), M(0), Nnz(0) {}
582 
583  /// Copy constructor
584  SparseMatrix(const SparseMatrix& source_matrix)
585  {
586  // Number of nonzero entries
587  Nnz = source_matrix.nnz();
588 
589  // Number of rows
590  N = source_matrix.nrow();
591 
592  // Number of columns
593  M = source_matrix.ncol();
594 
595  // Values stored in C-style array
596  Value = new T[Nnz];
597 
598  // Assign the values
599  for (unsigned long i = 0; i < Nnz; i++)
600  {
601  Value[i] = source_matrix.value()[i];
602  }
603  }
604 
605  /// Broken assignment operator
606  void operator=(const SparseMatrix&) = delete;
607 
608  /// Destructor, delete the memory associated with the values
609  virtual ~SparseMatrix()
610  {
611  delete[] Value;
612  Value = 0;
613  }
614 
615  /// Access to C-style value array
616  T* value()
617  {
618  return Value;
619  }
620 
621  /// Access to C-style value array (const version)
622  const T* value() const
623  {
624  return Value;
625  }
626 
627  /// Return the number of rows of the matrix
628  inline unsigned long nrow() const
629  {
630  return N;
631  }
632 
633  /// Return the number of columns of the matrix
634  inline unsigned long ncol() const
635  {
636  return M;
637  }
638 
639  /// Return the number of nonzero entries
640  inline unsigned long nnz() const
641  {
642  return Nnz;
643  }
644 
645  /// Output the "bottom right" entry regardless of it being
646  /// zero or not (this allows automatic detection of matrix size in
647  /// e.g. matlab, python).
648  virtual void output_bottom_right_zero_helper(std::ostream& outfile) const
649  {
650  std::string error_message = "SparseMatrix::output_bottom_right_zero_"
651  "helper() is a virtual function.\n";
652  error_message +=
653  "It must be overloaded for specific sparse matrix storage formats\n";
654 
655  throw OomphLibError(
656  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
657  }
658 
659  /// Indexed output function to print a matrix to the
660  /// stream outfile as i,j,a(i,j) for a(i,j)!=0 only.
661  virtual void sparse_indexed_output_helper(std::ostream& outfile) const
662  {
663  std::string error_message =
664  "SparseMatrix::sparse_indexed_output_helper() is a virtual function.\n";
665  error_message +=
666  "It must be overloaded for specific sparse matrix storage formats\n";
667 
668  throw OomphLibError(
669  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
670  }
671  };
672 
673 
674  //======================================================================
675  /// A class for compressed row matrices, a sparse storage format
676  /// Once again the recursive template trick is used to inform that base
677  /// class that is should use the access functions provided in the
678  /// CRMatrix class.
679  //=====================================================================
680  template<class T>
681  class CRMatrix : public SparseMatrix<T, CRMatrix<T>>
682  {
683  public:
684  /// Default constructor
686  {
687  Column_index = 0;
688  Row_start = 0;
689  }
690 
691 
692  /// Constructor: Pass vector of values, vector of column indices,
693  /// vector of row starts and number of rows and columns
694  /// Number of nonzero entries is read
695  /// off from value, so make sure the vector has been shrunk
696  /// to its correct length.
698  const Vector<int>& column_index_,
699  const Vector<int>& row_start_,
700  const unsigned long& n,
701  const unsigned long& m)
702  : SparseMatrix<T, CRMatrix<T>>()
703  {
704  Column_index = 0;
705  Row_start = 0;
706  build(value, column_index_, row_start_, n, m);
707  }
708 
709  /// Copy constructor
710  CRMatrix(const CRMatrix& source_matrix)
711  : SparseMatrix<T, CRMatrix<T>>(source_matrix)
712  {
713  // NNz, N and M are set the the copy constructor of the SparseMatrix
714  // called above
715  // Column indices stored in C-style array
716  Column_index = new int[this->Nnz];
717 
718  // Assign:
719  for (unsigned long i = 0; i < this->Nnz; i++)
720  {
721  Column_index[i] = source_matrix.column_index()[i];
722  }
723 
724  // Row start:
725  Row_start = new int[this->N + 1];
726 
727  // Assign:
728  for (unsigned long i = 0; i <= this->N; i++)
729  {
730  Row_start[i] = source_matrix.row_start()[i];
731  }
732  }
733 
734  /// Broken assignment operator
735  void operator=(const CRMatrix&) = delete;
736 
737  /// Destructor, delete any allocated memory
738  virtual ~CRMatrix()
739  {
740  delete[] Column_index;
741  Column_index = 0;
742  delete[] Row_start;
743  Row_start = 0;
744  }
745 
746  /// Access function that will be called by the read-only
747  /// round-bracket operator (const)
748  T get_entry(const unsigned long& i, const unsigned long& j) const
749  {
750 #ifdef RANGE_CHECKING
751  this->range_check(i, j);
752 #endif
753  for (long k = Row_start[i]; k < Row_start[i + 1]; k++)
754  {
755  if (unsigned(Column_index[k]) == j)
756  {
757  return this->Value[k];
758  }
759  }
760  return this->Zero;
761  }
762 
763  /// The read-write access function is deliberately broken
764  T& entry(const unsigned long& i, const unsigned long& j)
765  {
766  std::string error_string =
767  "Non-const access not provided for the CRMatrix<T> class\n";
768  error_string +=
769  "It is not possible to use round-bracket access: M(i,j)\n";
770  error_string += "if M is not declared as const.\n";
771  error_string += "The solution (albeit ugly) is to create const reference "
772  "to the matrix\n";
773  error_string += " const CRMatrix<T>& read_M = M;\n";
774  error_string += "Then read_M(i,j) is permitted\n";
775 
776  throw OomphLibError(
777  error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
778 
779  // Dummy return
780  T dummy;
781  return dummy;
782  }
783 
784  /// Access to C-style row_start array
785  int* row_start()
786  {
787  return Row_start;
788  }
789 
790  /// Access to C-style row_start array (const version)
791  const int* row_start() const
792  {
793  return Row_start;
794  }
795 
796  /// Access to C-style column index array
798  {
799  return Column_index;
800  }
801 
802  /// Access to C-style column index array (const version)
803  const int* column_index() const
804  {
805  return Column_index;
806  }
807 
808  /// Output the "bottom right" entry regardless of it being
809  /// zero or not (this allows automatic detection of matrix size in
810  /// e.g. matlab, python).
811  void output_bottom_right_zero_helper(std::ostream& outfile) const
812  {
813  int last_row_local = this->N - 1;
814  int last_col = this->M - 1;
815 
816  // Use this strange thingy because of the CRTP discussed above.
817  T last_value = this->operator()(last_row_local, last_col);
818 
819  if (last_value == T(0))
820  {
821  outfile << last_row_local << " " << last_col << " " << T(0)
822  << std::endl;
823  }
824  }
825 
826  /// Indexed output function to print a matrix to the
827  /// stream outfile as i,j,a(i,j) for a(i,j)!=0 only.
828  void sparse_indexed_output_helper(std::ostream& outfile) const
829  {
830  for (unsigned long i = 0; i < this->N; i++)
831  {
832  for (long j = Row_start[i]; j < Row_start[i + 1]; j++)
833  {
834  outfile << i << " " << Column_index[j] << " " << this->Value[j]
835  << std::endl;
836  }
837  }
838  }
839 
840  /// Wipe matrix data and set all values to 0.
842 
843  /// Build matrix from compressed representation. Number of nonzero
844  /// entries is read off from value, so make sure the vector has been shrunk
845  /// to its correct length. This matrix forms the storage for
846  /// CRDoubleMatrices which are distributable. The argument n should be the
847  /// number of local rows. The argument m is the number of columns
848  void build(const Vector<T>& value,
849  const Vector<int>& column_index,
850  const Vector<int>& row_start,
851  const unsigned long& n,
852  const unsigned long& m);
853 
854 
855  /// Function to build matrix from pointers to arrays
856  /// which hold the row starts, column indices and non-zero values.
857  /// The final two arguments are the number of rows and columns.
858  /// Note that, as the name suggests, this function does not
859  /// make a copy of the data pointed to by the first three arguments!
861  int* column_index,
862  int* row_start,
863  const unsigned long& nnz,
864  const unsigned long& n,
865  const unsigned long& m);
866 
867 
868  protected:
869  /// Column index
871 
872  /// Start index for row
873  int* Row_start;
874  };
875 
876 
877  // Forward definition for the superlu solver
878  class SuperLUSolver;
879 
880 
881  //=============================================================================
882  /// A class for compressed row matrices. This is a distributable
883  /// object.
884  //=============================================================================
885  class CRDoubleMatrix : public Matrix<double, CRDoubleMatrix>,
886  public DoubleMatrixBase,
888  {
889  public:
890  /// Default constructor
891  CRDoubleMatrix();
892 
893  /// Constructor: vector of values, vector of column indices,
894  /// vector of row starts and number of rows and columns.
896  const unsigned& ncol,
897  const Vector<double>& value,
898  const Vector<int>& column_index,
899  const Vector<int>& row_start);
900 
901  /// Constructor: just stores the distribution but does not build the
902  /// matrix
904 
905  /// Copy constructor
906  CRDoubleMatrix(const CRDoubleMatrix& matrix);
907 
908  /// Broken assignment operator
909  void operator=(const CRDoubleMatrix&) = delete;
910 
911  /// Destructor
912  virtual ~CRDoubleMatrix();
913 
914  /// Access function: returns the vector Index_of_diagonal_entries.
915  /// The i-th entry of the vector contains the index of the last entry
916  /// below or on the diagonal. If there are no entries below or on the
917  /// diagonal then the corresponding entry is -1. If, however, there are
918  /// no entries in the row then the entry is irrelevant and is kept
919  /// as the initialised value; 0.
921  {
922  // Check to see if the vector has been set up
923  if (Index_of_diagonal_entries.size() == 0)
924  {
925  // Make the warning
926  std::string err_strng =
927  "The Index_of_diagonal_entries vector has not been ";
928  err_strng += "set up yet. Run sort_entries() to set this vector up.";
929 
930  // Throw the warning
931  OomphLibWarning(err_strng,
932  "CRDoubleMatrix::get_index_of_diagonal_entries()",
933  OOMPH_EXCEPTION_LOCATION);
934  }
935 
936  // Return the vector
938  } // End of index_of_diagonal_entries
939 
940  /// Create a struct to provide a comparison function for std::sort
942  {
943  // Define the comparison operator
944  bool operator()(const std::pair<int, double>& pair_1,
945  const std::pair<int, double>& pair_2)
946  {
947  // If the first argument of pair_1 is less than the first argument of
948  // pair_2 then return TRUE otherwise return FALSE
949  return (pair_1.first < pair_2.first);
950  }
952 
953  /// Runs through the column index vector and checks if the entries
954  /// follow the regular lexicographical ordering of matrix entries, i.e.
955  /// it will check (at the i-th row of the matrix) if the entries in the
956  /// column index vector associated with this row are in increasing order
957  bool entries_are_sorted(const bool& doc_unordered_entries = false) const;
958 
959  /// Sorts the entries associated with each row of the matrix in the
960  /// column index vector and the value vector into ascending order and sets
961  /// up the Index_of_diagonal_entries vector
962  void sort_entries();
963 
964  /// build method: vector of values, vector of column indices,
965  /// vector of row starts and number of rows and columns.
967  const unsigned& ncol,
968  const Vector<double>& value,
969  const Vector<int>& column_index,
970  const Vector<int>& row_start);
971 
972  /// rebuild the matrix - assembles an empty matrix will a defined
973  /// distribution
975 
976  /// keeps the existing distribution and just matrix that is stored
977  void build(const unsigned& ncol,
978  const Vector<double>& value,
979  const Vector<int>& column_index,
980  const Vector<int>& row_start);
981 
982  /// keeps the existing distribution and just matrix that is stored
983  /// without copying the matrix data
984  void build_without_copy(const unsigned& ncol,
985  const unsigned& nnz,
986  double* value,
987  int* column_index,
988  int* row_start);
989 
990  /// The contents of the matrix are redistributed to match the new
991  /// distribution. In a non-MPI build this method does nothing.
992  /// \b NOTE 1: The current distribution and the new distribution must have
993  /// the same number of global rows.
994  /// \b NOTE 2: The current distribution and the new distribution must have
995  /// the same Communicator.
996  void redistribute(const LinearAlgebraDistribution* const& dist_pt);
997 
998  /// clear
999  void clear();
1000 
1001  /// Return the number of rows of the matrix
1002  inline unsigned long nrow() const
1003  {
1005  }
1006 
1007  /// Return the number of columns of the matrix
1008  inline unsigned long ncol() const
1009  {
1010  return CR_matrix.ncol();
1011  }
1012 
1013  /// Output the "bottom right" entry regardless of it being
1014  /// zero or not (this allows automatic detection of matrix size in
1015  /// e.g. matlab, python).
1016  void output_bottom_right_zero_helper(std::ostream& outfile) const
1017  {
1019  }
1020 
1021  /// Indexed output function to print a matrix to the
1022  /// stream outfile as i,j,a(i,j) for a(i,j)!=0 only.
1023  void sparse_indexed_output_helper(std::ostream& outfile) const
1024  {
1026  }
1027 
1028  /// Indexed output function to print a matrix to a
1029  /// file as i,j,a(i,j) for a(i,j)!=0 only. Specify filename.
1030  /// This uses acual global row numbers.
1032  {
1033  // Get offset
1034  unsigned first_row = distribution_pt()->first_row();
1035 
1036  // Open file
1037  std::ofstream some_file;
1038  some_file.open(filename.c_str());
1039  unsigned n = nrow_local();
1040  for (unsigned long i = 0; i < n; i++)
1041  {
1042  for (long j = row_start()[i]; j < row_start()[i + 1]; j++)
1043  {
1044  some_file << first_row + i << " " << column_index()[j] << " "
1045  << value()[j] << std::endl;
1046  }
1047  }
1048  some_file.close();
1049  }
1050 
1051  /// Overload the round-bracket access operator for read-only access. In a
1052  /// distributed matrix i refers to the local row index.
1053  inline double operator()(const unsigned long& i,
1054  const unsigned long& j) const
1055  {
1056  return CR_matrix.get_entry(i, j);
1057  }
1058 
1059  /// Access to C-style row_start array
1060  int* row_start()
1061  {
1062  return CR_matrix.row_start();
1063  }
1064 
1065  /// Access to C-style row_start array (const version)
1066  const int* row_start() const
1067  {
1068  return CR_matrix.row_start();
1069  }
1070 
1071  /// Access to C-style column index array
1073  {
1074  return CR_matrix.column_index();
1075  }
1076 
1077  /// Access to C-style column index array (const version)
1078  const int* column_index() const
1079  {
1080  return CR_matrix.column_index();
1081  }
1082 
1083  /// Access to C-style value array
1084  double* value()
1085  {
1086  return CR_matrix.value();
1087  }
1088 
1089  /// Access to C-style value array (const version)
1090  const double* value() const
1091  {
1092  return CR_matrix.value();
1093  }
1094 
1095  /// Return the number of nonzero entries (the local nnz)
1096  inline unsigned long nnz() const
1097  {
1098  return CR_matrix.nnz();
1099  }
1100 
1101  /// LU decomposition using SuperLU if matrix is not distributed or
1102  /// distributed onto a single processor.
1103  virtual void ludecompose();
1104 
1105  /// LU back solve for given RHS
1106  virtual void lubksub(DoubleVector& rhs);
1107 
1108  /// Multiply the matrix by the vector x: soln=Ax
1109  void multiply(const DoubleVector& x, DoubleVector& soln) const;
1110 
1111  /// Multiply the transposed matrix by the vector x: soln=A^T x
1112  void multiply_transpose(const DoubleVector& x, DoubleVector& soln) const;
1113 
1114  /// Function to multiply this matrix by the CRDoubleMatrix matrix_in.
1115  /// In a serial matrix, there are 4 methods available:
1116  /// Method 1: First runs through this matrix and matrix_in to find the
1117  /// storage
1118  /// requirements for result - arrays of the correct size are
1119  /// then allocated before performing the calculation.
1120  /// Minimises memory requirements but more costly.
1121  /// Method 2: Grows storage for values and column indices of result 'on the
1122  /// fly' using an array of maps. Faster but more memory
1123  /// intensive.
1124  /// Method 3: Grows storage for values and column indices of result 'on the
1125  /// fly' using a vector of vectors. Not particularly impressive
1126  /// on the platforms we tried...
1127  /// Method 4: Trilinos Epetra Matrix Matrix multiply.
1128  /// Method 5: Trilinox Epetra Matrix Matrix Mulitply (ml based)
1129  /// If Trilinos is installed then Method 4 is employed by default, otherwise
1130  /// Method 2 is employed by default.
1131  /// In a distributed matrix, only Trilinos Epetra Matrix Matrix multiply
1132  /// is available.
1133  void multiply(const CRDoubleMatrix& matrix_in,
1134  CRDoubleMatrix& result) const;
1135 
1136  /// For every row, find the maximum absolute value of the
1137  /// entries in this row. Set all values that are less than alpha times
1138  /// this maximum to zero and return the resulting matrix in
1139  /// reduced_matrix. Note: Diagonal entries are retained regardless
1140  /// of their size.
1141  void matrix_reduction(const double& alpha, CRDoubleMatrix& reduced_matrix);
1142 
1143  /// Access function to Serial_matrix_matrix_multiply_method, the flag
1144  /// which determines the matrix matrix multiplication method used for serial
1145  /// matrices.
1146  /// Method 1: First runs through this matrix and matrix_in to find the
1147  /// storage
1148  /// requirements for result - arrays of the correct size are
1149  /// then allocated before performing the calculation.
1150  /// Minimises memory requirements but more costly.
1151  /// Method 2: Grows storage for values and column indices of result 'on the
1152  /// fly' using an array of maps. Faster but more memory
1153  /// intensive.
1154  /// Method 3: Grows storage for values and column indices of result 'on the
1155  /// fly' using a vector of vectors. Not particularly impressive
1156  /// on the platforms we tried...
1157  /// Method 4: Trilinos Epetra Matrix Matrix multiply.
1158  /// Method 5: Trilinos Epetra Matrix Matrix multiply (ML based).
1160  {
1162  }
1163 
1164  /// Read only access function (const version) to
1165  /// Serial_matrix_matrix_multiply_method, the flag
1166  /// which determines the matrix matrix multiplication method used for serial
1167  /// matrices.
1168  /// Method 1: First runs through this matrix and matrix_in to find the
1169  /// storage
1170  /// requirements for result - arrays of the correct size are
1171  /// then allocated before performing the calculation.
1172  /// Minimises memory requirements but more costly.
1173  /// Method 2: Grows storage for values and column indices of result 'on the
1174  /// fly' using an array of maps. Faster but more memory
1175  /// intensive.
1176  /// Method 3: Grows storage for values and column indices of result 'on the
1177  /// fly' using a vector of vectors. Not particularly impressive
1178  /// on the platforms we tried...
1179  /// Method 4: Trilinos Epetra Matrix Matrix multiply.
1180  /// Method 5: Trilinos Epetra Matrix Matrix multiply (ML based).
1182  {
1184  }
1185 
1186  /// Access function to Distributed_matrix_matrix_multiply_method, the
1187  /// flag which determines the matrix matrix multiplication method used for
1188  /// distributed matrices.
1189  /// Method 1: Trilinos Epetra Matrix Matrix multiply.
1190  /// Method 2: Trilinos Epetra Matrix Matrix multiply (ML based).
1192  {
1194  }
1195 
1196  /// Read only access function (const version) to
1197  /// Distributed_matrix_matrix_multiply_method, the
1198  /// flag which determines the matrix matrix multiplication method used for
1199  /// distributed matrices.
1200  /// Method 1: Trilinos Epetra Matrix Matrix multiply.
1201  /// Method 2: Trilinos Epetra Matrix Matrix multiply (ML based).
1203  {
1205  }
1206 
1207  /// access function to the Built flag - indicates whether the matrix
1208  /// has been build - i.e. the distribution has been defined and the matrix
1209  /// assembled.
1210  bool built() const
1211  {
1212  return Built;
1213  }
1214 
1215  /// if this matrix is distributed then a the equivalent global matrix
1216  /// is built using new and returned. The calling method is responsible for
1217  /// the destruction of the new matrix.
1218  CRDoubleMatrix* global_matrix() const;
1219 
1220  /// Returns the transpose of this matrix
1221  void get_matrix_transpose(CRDoubleMatrix* result) const;
1222 
1223  /// returns the inf-norm of this matrix
1224  double inf_norm() const;
1225 
1226  /// returns a Vector of diagonal entries of this matrix.
1227  /// This only works with square matrices. This condition may be relaxed
1228  /// in the future if need be.
1230 
1231  /// element-wise addition of this matrix with matrix_in.
1232  void add(const CRDoubleMatrix& matrix_in,
1233  CRDoubleMatrix& result_matrix) const;
1234 
1235  private:
1236  /// Vector whose i'th entry contains the index of the last entry
1237  /// below or on the diagonal of the i'th row of the matrix
1239 
1240  /// Flag to determine which matrix-matrix multiplication method is
1241  /// used (for serial (or global) matrices)
1243 
1244  /// Flag to determine which matrix-matrix multiplication method is
1245  /// used (for distributed matrices)
1247 
1248  /// Storage for the Matrix in CR Format
1250 
1251  /// Flag to indicate whether the matrix has been built - i.e. the
1252  /// distribution has been setup AND the matrix has been assembled.
1253  bool Built;
1254  };
1255 
1256 
1257  /// ////////////////////////////////////////////////////////////////////////////
1258  /// ////////////////////////////////////////////////////////////////////////////
1259  /// ////////////////////////////////////////////////////////////////////////////
1260 
1261 
1262  // Forward definition of the DenseLU class
1263  class DenseLU;
1264 
1265  //=================================================================
1266  /// Class of matrices containing doubles, and stored as a
1267  /// DenseMatrix<double>, but with solving functionality inherited
1268  /// from the abstract DoubleMatrix class.
1269  //=================================================================
1270  class DenseDoubleMatrix : public DoubleMatrixBase, public DenseMatrix<double>
1271  {
1272  public:
1273  /// Constructor, set the default linear solver
1275 
1276  /// Constructor to build a square n by n matrix.
1277  DenseDoubleMatrix(const unsigned long& n);
1278 
1279  /// Constructor to build a matrix with n rows and m columns.
1280  DenseDoubleMatrix(const unsigned long& n, const unsigned long& m);
1281 
1282  /// Constructor to build a matrix with n rows and m columns,
1283  /// with initial value initial_val
1284  DenseDoubleMatrix(const unsigned long& n,
1285  const unsigned long& m,
1286  const double& initial_val);
1287 
1288  /// Broken copy constructor
1289  DenseDoubleMatrix(const DenseDoubleMatrix& matrix) = delete;
1290 
1291  /// Broken assignment operator
1292  void operator=(const DenseDoubleMatrix&) = delete;
1293 
1294  /// Return the number of rows of the matrix
1295  inline unsigned long nrow() const
1296  {
1297  return DenseMatrix<double>::nrow();
1298  }
1299 
1300  /// Return the number of columns of the matrix
1301  inline unsigned long ncol() const
1302  {
1303  return DenseMatrix<double>::ncol();
1304  }
1305 
1306  /// Overload the const version of the round-bracket access operator
1307  /// for read-only access.
1308  inline double operator()(const unsigned long& i,
1309  const unsigned long& j) const
1310  {
1311  return DenseMatrix<double>::get_entry(i, j);
1312  }
1313 
1314  /// Overload the non-const version of the round-bracket access
1315  /// operator for read-write access
1316  inline double& operator()(const unsigned long& i, const unsigned long& j)
1317  {
1318  return DenseMatrix<double>::entry(i, j);
1319  }
1320 
1321  /// Destructor
1322  virtual ~DenseDoubleMatrix();
1323 
1324  /// LU decomposition using DenseLU (default linea solver)
1325  virtual void ludecompose();
1326 
1327  /// LU backsubstitution
1328  virtual void lubksub(DoubleVector& rhs);
1329 
1330  /// LU backsubstitution
1331  virtual void lubksub(Vector<double>& rhs);
1332 
1333  /// Determine eigenvalues and eigenvectors, using
1334  /// Jacobi rotations. Only for symmetric matrices. Nothing gets overwritten!
1335  /// - \c eigen_vect(i,j) = j-th component of i-th eigenvector.
1336  /// - \c eigen_val(i) is the i-th eigenvalue; same ordering as in
1337  /// eigenvectors
1338  void eigenvalues_by_jacobi(Vector<double>& eigen_val,
1339  DenseMatrix<double>& eigen_vect) const;
1340 
1341  /// Multiply the matrix by the vector x: soln=Ax
1342  void multiply(const DoubleVector& x, DoubleVector& soln) const;
1343 
1344  /// Multiply the transposed matrix by the vector x: soln=A^T x
1345  void multiply_transpose(const DoubleVector& x, DoubleVector& soln) const;
1346 
1347  /// For every row, find the maximum absolute value of the
1348  /// entries in this row. Set all values that are less than alpha times
1349  /// this maximum to zero and return the resulting matrix in
1350  /// reduced_matrix. Note: Diagonal entries are retained regardless
1351  /// of their size.
1352  void matrix_reduction(const double& alpha,
1353  DenseDoubleMatrix& reduced_matrix);
1354 
1355  /// Function to multiply this matrix by a DenseDoubleMatrix matrix_in
1356  void multiply(const DenseDoubleMatrix& matrix_in,
1357  DenseDoubleMatrix& result);
1358  };
1359 
1360  /// //////////////////////////////////////////////////////////////////
1361  /// //////////////////////////////////////////////////////////////////
1362  /// //////////////////////////////////////////////////////////////////
1363 
1364 
1365  //=================================================================
1366  /// A Rank 3 Tensor class
1367  //=================================================================
1368  template<class T>
1370  {
1371  private:
1372  /// Private internal representation as pointer to data
1374 
1375  /// 1st Tensor dimension
1376  unsigned N;
1377 
1378  /// 2nd Tensor dimension
1379  unsigned M;
1380 
1381  /// 3rd Tensor dimension
1382  unsigned P;
1383 
1384  /// Range check to catch when an index is out of bounds, if so, it
1385  /// issues a warning message and dies by throwing an \c OomphLibError
1386  void range_check(const unsigned long& i,
1387  const unsigned long& j,
1388  const unsigned long& k) const
1389  {
1390  if (i >= N)
1391  {
1392  std::ostringstream error_message;
1393  error_message << "Range Error: i=" << i << " is not in the range (0,"
1394  << N - 1 << ")." << std::endl;
1395 
1396  throw OomphLibError(error_message.str(),
1397  OOMPH_CURRENT_FUNCTION,
1398  OOMPH_EXCEPTION_LOCATION);
1399  }
1400  else if (j >= M)
1401  {
1402  std::ostringstream error_message;
1403  error_message << "Range Error: j=" << j << " is not in the range (0,"
1404  << M - 1 << ")." << std::endl;
1405 
1406  throw OomphLibError(error_message.str(),
1407  OOMPH_CURRENT_FUNCTION,
1408  OOMPH_EXCEPTION_LOCATION);
1409  }
1410  else if (k >= P)
1411  {
1412  std::ostringstream error_message;
1413  error_message << "Range Error: k=" << k << " is not in the range (0,"
1414  << P - 1 << ")." << std::endl;
1415 
1416  throw OomphLibError(error_message.str(),
1417  OOMPH_CURRENT_FUNCTION,
1418  OOMPH_EXCEPTION_LOCATION);
1419  }
1420  }
1421 
1422 
1423  public:
1424  /// Empty constructor
1425  RankThreeTensor() : Tensordata(0), N(0), M(0), P(0) {}
1426 
1427  /// Copy constructor: Deep copy
1428  RankThreeTensor(const RankThreeTensor& source_tensor)
1429  {
1430  // Set row and column lengths
1431  N = source_tensor.nindex1();
1432  M = source_tensor.nindex2();
1433  P = source_tensor.nindex3();
1434  // Assign space for the data
1435  Tensordata = new T[N * M * P];
1436  // Copy the data across from the other matrix
1437  for (unsigned i = 0; i < N; i++)
1438  {
1439  for (unsigned j = 0; j < M; j++)
1440  {
1441  for (unsigned k = 0; k < P; k++)
1442  {
1443  Tensordata[P * (M * i + j) + k] = source_tensor(i, j, k);
1444  }
1445  }
1446  }
1447  }
1448 
1449  /// Copy assignement
1451  {
1452  // Don't create a new matrix if the assignement is the identity
1453  if (this != &source_tensor)
1454  {
1455  // Check row and column length
1456  unsigned long n = source_tensor.nindex1();
1457  unsigned long m = source_tensor.nindex2();
1458  unsigned long p = source_tensor.nindex3();
1459  // Resie the tensor to be the same size as the old tensor
1460  if ((N != n) || (M != m) || (P != p))
1461  {
1462  resize(n, m, p);
1463  }
1464 
1465  // Copy entries across from the other matrix
1466  for (unsigned long i = 0; i < N; i++)
1467  {
1468  for (unsigned long j = 0; j < M; j++)
1469  {
1470  for (unsigned long k = 0; k < P; k++)
1471  {
1472  (*this)(i, j, k) = source_tensor(i, j, k);
1473  }
1474  }
1475  }
1476  }
1477  // Return reference to object itself (i.e. de-reference this pointer)
1478  return *this;
1479  }
1480 
1481 
1482  /// One parameter constructor produces a cubic nxnxn tensor
1483  RankThreeTensor(const unsigned long& n)
1484  {
1485  // Set row and column lengths
1486  N = n;
1487  M = n;
1488  P = n;
1489  // Assign space for the n rows
1490  Tensordata = new T[N * M * P];
1491  // Initialise to zero if required
1492 #ifdef OOMPH_INITIALISE_DENSE_MATRICES
1493  initialise(T(0));
1494 #endif
1495  }
1496 
1497  /// Three parameter constructor, general non-square tensor
1498  RankThreeTensor(const unsigned long& n_index1,
1499  const unsigned long& n_index2,
1500  const unsigned long& n_index3)
1501  {
1502  // Set row and column lengths
1503  N = n_index1;
1504  M = n_index2;
1505  P = n_index3;
1506  // Assign space for the n rows
1507  Tensordata = new T[N * M * P];
1508  // Initialise to zero if required
1509 #ifdef OOMPH_INITIALISE_DENSE_MATRICES
1510  initialise(T(0));
1511 #endif
1512  }
1513 
1514 
1515  /// Three parameter constructor, general non-square tensor
1516  RankThreeTensor(const unsigned long& n_index1,
1517  const unsigned long& n_index2,
1518  const unsigned long& n_index3,
1519  const T& initial_val)
1520  {
1521  // Set row and column lengths
1522  N = n_index1;
1523  M = n_index2;
1524  P = n_index3;
1525  // Assign space for the n rows
1526  Tensordata = new T[N * M * P];
1527  // Initialise to the initial value
1528  initialise(initial_val);
1529  }
1530 
1531  /// Destructor: delete the pointers
1533  {
1534  delete[] Tensordata;
1535  Tensordata = 0;
1536  }
1537 
1538  /// Resize to a square nxnxn tensor
1539  void resize(const unsigned long& n)
1540  {
1541  resize(n, n, n);
1542  }
1543 
1544  /// Resize to a general tensor
1545  void resize(const unsigned long& n_index1,
1546  const unsigned long& n_index2,
1547  const unsigned long& n_index3)
1548  {
1549  // If the sizes have not changed do nothing
1550  if ((n_index1 == N) && (n_index2 == M) && (n_index3 == P))
1551  {
1552  return;
1553  }
1554  // Store old sizes
1555  unsigned long n_old = N, m_old = M, p_old = P;
1556  // Reassign the sizes
1557  N = n_index1;
1558  M = n_index2;
1559  P = n_index3;
1560  // Store triple pointer to old matrix data
1561  T* temp_tensor = Tensordata;
1562  // Re-create Tensordata in new size
1563  Tensordata = new T[N * M * P];
1564 #ifdef OOMPH_INITIALISE_DENSE_MATRICES
1565  initialise(T(0));
1566 #endif
1567  // Transfer values
1568  unsigned long n_copy, m_copy, p_copy;
1569  n_copy = std::min(n_old, n_index1);
1570  m_copy = std::min(m_old, n_index2);
1571  p_copy = std::min(p_old, n_index3);
1572  // If matrix has values, transfer them to new matrix
1573  // Loop over rows
1574  for (unsigned long i = 0; i < n_copy; i++)
1575  {
1576  // Loop over columns
1577  for (unsigned long j = 0; j < m_copy; j++)
1578  {
1579  // Loop over columns
1580  for (unsigned long k = 0; k < p_copy; k++)
1581  {
1582  // Transfer values from temp_tensor
1583  Tensordata[M * P * i + P * j + k] =
1584  temp_tensor[m_old * p_old * i + p_old * j + k];
1585  }
1586  }
1587  }
1588  // Now kill storage for old tensor
1589  delete[] temp_tensor;
1590  }
1591 
1592  /// Resize to a general tensor
1593  void resize(const unsigned long& n_index1,
1594  const unsigned long& n_index2,
1595  const unsigned long& n_index3,
1596  const T& initial_value)
1597  {
1598  // If the sizes have not changed do nothing
1599  if ((n_index1 == N) && (n_index2 == M) && (n_index3 == P))
1600  {
1601  return;
1602  }
1603  // Store old sizes
1604  unsigned long n_old = N, m_old = M, p_old = P;
1605  // Reassign the sizes
1606  N = n_index1;
1607  M = n_index2;
1608  P = n_index3;
1609  // Store triple pointer to old matrix data
1610  T* temp_tensor = Tensordata;
1611  // Re-create Tensordata in new size
1612  Tensordata = new T[N * M * P];
1613  // Initialise the newly allocated storage
1614  initialise(initial_value);
1615 
1616  // Transfer values
1617  unsigned long n_copy, m_copy, p_copy;
1618  n_copy = std::min(n_old, n_index1);
1619  m_copy = std::min(m_old, n_index2);
1620  p_copy = std::min(p_old, n_index3);
1621  // If matrix has values, transfer them to new matrix
1622  // Loop over rows
1623  for (unsigned long i = 0; i < n_copy; i++)
1624  {
1625  // Loop over columns
1626  for (unsigned long j = 0; j < m_copy; j++)
1627  {
1628  // Loop over columns
1629  for (unsigned long k = 0; k < p_copy; k++)
1630  {
1631  // Transfer values from temp_tensor
1632  Tensordata[M * P * i + P * j + k] =
1633  temp_tensor[m_old * p_old * i + p_old * j + k];
1634  }
1635  }
1636  }
1637  // Now kill storage for old tensor
1638  delete[] temp_tensor;
1639  }
1640 
1641  /// Initialise all values in the tensor to val
1642  void initialise(const T& val)
1643  {
1644  for (unsigned long i = 0; i < (N * M * P); ++i)
1645  {
1646  Tensordata[i] = val;
1647  }
1648  }
1649 
1650  /// Return the range of index 1 of the tensor
1651  unsigned long nindex1() const
1652  {
1653  return N;
1654  }
1655 
1656  /// Return the range of index 2 of the tensor
1657  unsigned long nindex2() const
1658  {
1659  return M;
1660  }
1661 
1662  /// Return the range of index 3 of the tensor
1663  unsigned long nindex3() const
1664  {
1665  return P;
1666  }
1667 
1668  /// Overload the round brackets to give access as a(i,j,k)
1669  inline T& operator()(const unsigned long& i,
1670  const unsigned long& j,
1671  const unsigned long& k)
1672  {
1673 #ifdef RANGE_CHECKING
1674  this->range_check(i, j, k);
1675 #endif
1676  return Tensordata[P * (M * i + j) + k];
1677  }
1678 
1679  /// Overload a const version for read-only access as a(i,j,k)
1680  inline T operator()(const unsigned long& i,
1681  const unsigned long& j,
1682  const unsigned long& k) const
1683  {
1684 #ifdef RANGE_CHECKING
1685  this->range_check(i, j, k);
1686 #endif
1687  return Tensordata[P * (M * i + j) + k];
1688  }
1689  };
1690 
1691  /// //////////////////////////////////////////////////////////////////
1692  /// //////////////////////////////////////////////////////////////////
1693  /// //////////////////////////////////////////////////////////////////
1694 
1695 
1696  //=================================================================
1697  /// A Rank 4 Tensor class
1698  //=================================================================
1699  template<class T>
1701  {
1702  private:
1703  /// Private internal representation as pointer to data
1705 
1706  /// 1st Tensor dimension
1707  unsigned N;
1708 
1709  /// 2nd Tensor dimension
1710  unsigned M;
1711 
1712  /// 3rd Tensor dimension
1713  unsigned P;
1714 
1715  /// 4th Tensor dimension
1716  unsigned Q;
1717 
1718  /// Range check to catch when an index is out of bounds, if so, it
1719  /// issues a warning message and dies by throwing an \c OomphLibError
1720  void range_check(const unsigned long& i,
1721  const unsigned long& j,
1722  const unsigned long& k,
1723  const unsigned long& l) const
1724  {
1725  if (i >= N)
1726  {
1727  std::ostringstream error_message;
1728  error_message << "Range Error: i=" << i << " is not in the range (0,"
1729  << N - 1 << ")." << std::endl;
1730 
1731  throw OomphLibError(error_message.str(),
1732  OOMPH_CURRENT_FUNCTION,
1733  OOMPH_EXCEPTION_LOCATION);
1734  }
1735  else if (j >= M)
1736  {
1737  std::ostringstream error_message;
1738  error_message << "Range Error: j=" << j << " is not in the range (0,"
1739  << M - 1 << ")." << std::endl;
1740 
1741  throw OomphLibError(error_message.str(),
1742  OOMPH_CURRENT_FUNCTION,
1743  OOMPH_EXCEPTION_LOCATION);
1744  }
1745  else if (k >= P)
1746  {
1747  std::ostringstream error_message;
1748  error_message << "Range Error: k=" << k << " is not in the range (0,"
1749  << P - 1 << ")." << std::endl;
1750 
1751  throw OomphLibError(error_message.str(),
1752  OOMPH_CURRENT_FUNCTION,
1753  OOMPH_EXCEPTION_LOCATION);
1754  }
1755  else if (l >= Q)
1756  {
1757  std::ostringstream error_message;
1758  error_message << "Range Error: l=" << l << " is not in the range (0,"
1759  << Q - 1 << ")." << std::endl;
1760 
1761  throw OomphLibError(error_message.str(),
1762  OOMPH_CURRENT_FUNCTION,
1763  OOMPH_EXCEPTION_LOCATION);
1764  }
1765  }
1766 
1767  public:
1768  /// Empty constructor
1769  RankFourTensor() : Tensordata(0), N(0), M(0), P(0), Q(0) {}
1770 
1771  /// Copy constructor: Deep copy
1772  RankFourTensor(const RankFourTensor& source_tensor)
1773  {
1774  // Set row and column lengths
1775  N = source_tensor.nindex1();
1776  M = source_tensor.nindex2();
1777  P = source_tensor.nindex3();
1778  Q = source_tensor.nindex4();
1779 
1780  // Assign space for the data
1781  Tensordata = new T[N * M * P * Q];
1782 
1783  // Copy the data across from the other matrix
1784  for (unsigned i = 0; i < N; i++)
1785  {
1786  for (unsigned j = 0; j < M; j++)
1787  {
1788  for (unsigned k = 0; k < P; k++)
1789  {
1790  for (unsigned l = 0; l < Q; l++)
1791  {
1792  Tensordata[Q * (P * (M * i + j) + k) + l] =
1793  source_tensor(i, j, k, l);
1794  }
1795  }
1796  }
1797  }
1798  }
1799 
1800  /// Copy assignement
1801  RankFourTensor& operator=(const RankFourTensor& source_tensor)
1802  {
1803  // Don't create a new matrix if the assignement is the identity
1804  if (this != &source_tensor)
1805  {
1806  // Check row and column length
1807  unsigned long n = source_tensor.nindex1();
1808  unsigned long m = source_tensor.nindex2();
1809  unsigned long p = source_tensor.nindex3();
1810  unsigned long q = source_tensor.nindex4();
1811  // Resize the tensor to be the same size as the old tensor
1812  if ((N != n) || (M != m) || (P != p) || (Q != q))
1813  {
1814  resize(n, m, p, q);
1815  }
1816 
1817  // Copy entries across from the other matrix
1818  for (unsigned long i = 0; i < N; i++)
1819  {
1820  for (unsigned long j = 0; j < M; j++)
1821  {
1822  for (unsigned long k = 0; k < P; k++)
1823  {
1824  for (unsigned long l = 0; l < Q; l++)
1825  {
1826  (*this)(i, j, k, l) = source_tensor(i, j, k, l);
1827  }
1828  }
1829  }
1830  }
1831  }
1832  // Return reference to object itself (i.e. de-reference this pointer)
1833  return *this;
1834  }
1835 
1836 
1837  /// One parameter constructor produces a nxnxnxn tensor
1838  RankFourTensor(const unsigned long& n)
1839  {
1840  // Set row and column lengths
1841  N = n;
1842  M = n;
1843  P = n;
1844  Q = n;
1845  // Assign space for the n rows
1846  Tensordata = new T[N * M * P * Q];
1847  // Initialise to zero if required
1848 #ifdef OOMPH_INITIALISE_DENSE_MATRICES
1849  initialise(T(0));
1850 #endif
1851  }
1852 
1853  /// Four parameter constructor, general non-square tensor
1854  RankFourTensor(const unsigned long& n_index1,
1855  const unsigned long& n_index2,
1856  const unsigned long& n_index3,
1857  const unsigned long& n_index4)
1858  {
1859  // Set row and column lengths
1860  N = n_index1;
1861  M = n_index2;
1862  P = n_index3;
1863  Q = n_index4;
1864  // Assign space for the n rows
1865  Tensordata = new T[N * M * P * Q];
1866  // Initialise to zero if required
1867 #ifdef OOMPH_INITIALISE_DENSE_MATRICES
1868  initialise(T(0));
1869 #endif
1870  }
1871 
1872 
1873  /// Four parameter constructor, general non-square tensor
1874  RankFourTensor(const unsigned long& n_index1,
1875  const unsigned long& n_index2,
1876  const unsigned long& n_index3,
1877  const unsigned long& n_index4,
1878  const T& initial_val)
1879  {
1880  // Set row and column lengths
1881  N = n_index1;
1882  M = n_index2;
1883  P = n_index3;
1884  Q = n_index4;
1885  // Assign space for the n rows
1886  Tensordata = new T[N * M * P * Q];
1887  // Initialise to the initial value
1888  initialise(initial_val);
1889  }
1890 
1891  /// Destructor: delete the pointers
1893  {
1894  delete[] Tensordata;
1895  Tensordata = 0;
1896  }
1897 
1898  /// Resize to a square nxnxnxn tensor
1899  void resize(const unsigned long& n)
1900  {
1901  resize(n, n, n, n);
1902  }
1903 
1904  /// Resize to a general tensor
1905  void resize(const unsigned long& n_index1,
1906  const unsigned long& n_index2,
1907  const unsigned long& n_index3,
1908  const unsigned long& n_index4)
1909  {
1910  // If the sizes have not changed do nothing
1911  if ((n_index1 == N) && (n_index2 == M) && (n_index3 == P) &&
1912  (n_index4 == Q))
1913  {
1914  return;
1915  }
1916  // Store old sizes
1917  unsigned long n_old = N, m_old = M, p_old = P, q_old = Q;
1918  // Reassign the sizes
1919  N = n_index1;
1920  M = n_index2;
1921  P = n_index3;
1922  Q = n_index4;
1923  // Store pointer to old matrix data
1924  T* temp_tensor = Tensordata;
1925  // Re-create Tensordata in new size
1926  Tensordata = new T[N * M * P * Q];
1927 #ifdef OOMPH_INITIALISE_DENSE_MATRICES
1928  initialise(T(0));
1929 #endif
1930  // Transfer values
1931  unsigned long n_copy, m_copy, p_copy, q_copy;
1932  n_copy = std::min(n_old, n_index1);
1933  m_copy = std::min(m_old, n_index2);
1934  p_copy = std::min(p_old, n_index3);
1935  q_copy = std::min(q_old, n_index4);
1936  // If matrix has values, transfer them to new matrix
1937  // Loop over rows
1938  for (unsigned long i = 0; i < n_copy; i++)
1939  {
1940  // Loop over columns
1941  for (unsigned long j = 0; j < m_copy; j++)
1942  {
1943  // Loop over columns
1944  for (unsigned long k = 0; k < p_copy; k++)
1945  {
1946  // Loop over columns
1947  for (unsigned long l = 0; l < q_copy; l++)
1948  {
1949  // Transfer values from temp_tensor
1950  Tensordata[Q * (M * P * i + P * j + k) + l] =
1951  temp_tensor[q_old * (m_old * p_old * i + p_old * j + k) + l];
1952  }
1953  }
1954  }
1955  }
1956  // Now kill storage for old tensor
1957  delete[] temp_tensor;
1958  }
1959 
1960  /// Resize to a general tensor
1961  void resize(const unsigned long& n_index1,
1962  const unsigned long& n_index2,
1963  const unsigned long& n_index3,
1964  const unsigned long& n_index4,
1965  const T& initial_value)
1966  {
1967  // If the sizes have not changed do nothing
1968  if ((n_index1 == N) && (n_index2 == M) && (n_index3 == P) &&
1969  (n_index4 == Q))
1970  {
1971  return;
1972  }
1973  // Store old sizes
1974  unsigned long n_old = N, m_old = M, p_old = P, q_old = Q;
1975  // Reassign the sizes
1976  N = n_index1;
1977  M = n_index2;
1978  P = n_index3;
1979  Q = n_index4;
1980  // Store triple pointer to old matrix data
1981  T* temp_tensor = Tensordata;
1982  // Re-create Tensordata in new size
1983  Tensordata = new T[N * M * P * Q];
1984  // Initialise the newly allocated storage
1985  initialise(initial_value);
1986 
1987  // Transfer values
1988  unsigned long n_copy, m_copy, p_copy, q_copy;
1989  n_copy = std::min(n_old, n_index1);
1990  m_copy = std::min(m_old, n_index2);
1991  p_copy = std::min(p_old, n_index3);
1992  q_copy = std::min(q_old, n_index4);
1993  // If matrix has values, transfer them to new matrix
1994  // Loop over rows
1995  for (unsigned long i = 0; i < n_copy; i++)
1996  {
1997  // Loop over columns
1998  for (unsigned long j = 0; j < m_copy; j++)
1999  {
2000  // Loop over columns
2001  for (unsigned long k = 0; k < p_copy; k++)
2002  {
2003  // Loop over columns
2004  for (unsigned long l = 0; l < q_copy; l++)
2005  {
2006  // Transfer values from temp_tensor
2007  Tensordata[Q * (M * P * i + P * j + k) + l] =
2008  temp_tensor[q_old * (m_old * p_old * i + p_old * j + k) + l];
2009  }
2010  }
2011  }
2012  }
2013  // Now kill storage for old tensor
2014  delete[] temp_tensor;
2015  }
2016 
2017  /// Initialise all values in the tensor to val
2018  void initialise(const T& val)
2019  {
2020  for (unsigned long i = 0; i < (N * M * P * Q); ++i)
2021  {
2022  Tensordata[i] = val;
2023  }
2024  }
2025 
2026  /// Return the range of index 1 of the tensor
2027  unsigned long nindex1() const
2028  {
2029  return N;
2030  }
2031 
2032  /// Return the range of index 2 of the tensor
2033  unsigned long nindex2() const
2034  {
2035  return M;
2036  }
2037 
2038  /// Return the range of index 3 of the tensor
2039  unsigned long nindex3() const
2040  {
2041  return P;
2042  }
2043 
2044  /// Return the range of index 4 of the tensor
2045  unsigned long nindex4() const
2046  {
2047  return Q;
2048  }
2049 
2050  /// Overload the round brackets to give access as a(i,j,k,l)
2051  inline T& operator()(const unsigned long& i,
2052  const unsigned long& j,
2053  const unsigned long& k,
2054  const unsigned long& l)
2055  {
2056 #ifdef RANGE_CHECKING
2057  this->range_check(i, j, k, l);
2058 #endif
2059  return Tensordata[Q * (P * (M * i + j) + k) + l];
2060  }
2061 
2062  /// Overload a const version for read-only access as a(i,j,k,l)
2063  inline T operator()(const unsigned long& i,
2064  const unsigned long& j,
2065  const unsigned long& k,
2066  const unsigned long& l) const
2067  {
2068 #ifdef RANGE_CHECKING
2069  this->range_check(i, j, k, l);
2070 #endif
2071  return Tensordata[Q * (P * (M * i + j) + k) + l];
2072  }
2073 
2074  /// Direct access to internal storage of data in flat-packed C-style
2075  /// column-major format. WARNING: Only for experienced users. Only
2076  /// use this if raw speed is of the essence, as in the solid mechanics
2077  /// problems.
2078  inline T& raw_direct_access(const unsigned long& i)
2079  {
2080  return Tensordata[i];
2081  }
2082 
2083  /// Direct access to internal storage of data in flat-packed C-style
2084  /// column-major format. WARNING: Only for experienced users. Only
2085  /// use this if raw speed is of the essence, as in the solid mechanics
2086  /// problems.
2087  inline const T& raw_direct_access(const unsigned long& i) const
2088  {
2089  return Tensordata[i];
2090  }
2091 
2092  /// Caculate the offset in flat-packed C-style, column-major format,
2093  /// required for a given i,j. WARNING: Only for experienced users. Only
2094  /// use this if raw speed is of the essence, as in the solid mechanics
2095  /// problems.
2096  unsigned offset(const unsigned long& i, const unsigned long& j) const
2097  {
2098  return (Q * (P * (M * i + j) + 0) + 0);
2099  }
2100  };
2101 
2102 
2103  /// ///////////////////////////////////////////////////////////////
2104  /// ///////////////////////////////////////////////////////////////
2105  /// ///////////////////////////////////////////////////////////////
2106 
2107 
2108  //=================================================================
2109  /// A Rank 5 Tensor class
2110  //=================================================================
2111  template<class T>
2113  {
2114  private:
2115  /// Private internal representation as pointer to data
2117 
2118  /// 1st Tensor dimension
2119  unsigned N;
2120 
2121  /// 2nd Tensor dimension
2122  unsigned M;
2123 
2124  /// 3rd Tensor dimension
2125  unsigned P;
2126 
2127  /// 4th Tensor dimension
2128  unsigned Q;
2129 
2130  /// 5th Tensor dimension
2131  unsigned R;
2132 
2133  /// Range check to catch when an index is out of bounds, if so, it
2134  /// issues a warning message and dies by throwing an \c OomphLibError
2135  void range_check(const unsigned long& i,
2136  const unsigned long& j,
2137  const unsigned long& k,
2138  const unsigned long& l,
2139  const unsigned long& m) const
2140  {
2141  if (i >= N)
2142  {
2143  std::ostringstream error_message;
2144  error_message << "Range Error: i=" << i << " is not in the range (0,"
2145  << N - 1 << ")." << std::endl;
2146 
2147  throw OomphLibError(error_message.str(),
2148  OOMPH_CURRENT_FUNCTION,
2149  OOMPH_EXCEPTION_LOCATION);
2150  }
2151  else if (j >= M)
2152  {
2153  std::ostringstream error_message;
2154  error_message << "Range Error: j=" << j << " is not in the range (0,"
2155  << M - 1 << ")." << std::endl;
2156 
2157  throw OomphLibError(error_message.str(),
2158  OOMPH_CURRENT_FUNCTION,
2159  OOMPH_EXCEPTION_LOCATION);
2160  }
2161  else if (k >= P)
2162  {
2163  std::ostringstream error_message;
2164  error_message << "Range Error: k=" << k << " is not in the range (0,"
2165  << P - 1 << ")." << std::endl;
2166 
2167  throw OomphLibError(error_message.str(),
2168  OOMPH_CURRENT_FUNCTION,
2169  OOMPH_EXCEPTION_LOCATION);
2170  }
2171  else if (l >= Q)
2172  {
2173  std::ostringstream error_message;
2174  error_message << "Range Error: l=" << l << " is not in the range (0,"
2175  << Q - 1 << ")." << std::endl;
2176 
2177  throw OomphLibError(error_message.str(),
2178  OOMPH_CURRENT_FUNCTION,
2179  OOMPH_EXCEPTION_LOCATION);
2180  }
2181  else if (m >= R)
2182  {
2183  std::ostringstream error_message;
2184  error_message << "Range Error: m=" << m << " is not in the range (0,"
2185  << R - 1 << ")." << std::endl;
2186 
2187  throw OomphLibError(error_message.str(),
2188  OOMPH_CURRENT_FUNCTION,
2189  OOMPH_EXCEPTION_LOCATION);
2190  }
2191  }
2192 
2193  public:
2194  /// Empty constructor
2195  RankFiveTensor() : Tensordata(0), N(0), M(0), P(0), Q(0), R(0) {}
2196 
2197  /// Copy constructor: Deep copy
2198  RankFiveTensor(const RankFiveTensor& source_tensor)
2199  {
2200  // Set row and column lengths
2201  N = source_tensor.nindex1();
2202  M = source_tensor.nindex2();
2203  P = source_tensor.nindex3();
2204  Q = source_tensor.nindex4();
2205  R = source_tensor.nindex5();
2206 
2207  // Assign space for the data
2208  Tensordata = new T[N * M * P * Q * R];
2209 
2210  // Copy the data across from the other matrix
2211  for (unsigned i = 0; i < N; i++)
2212  {
2213  for (unsigned j = 0; j < M; j++)
2214  {
2215  for (unsigned k = 0; k < P; k++)
2216  {
2217  for (unsigned l = 0; l < Q; l++)
2218  {
2219  for (unsigned m = 0; m < R; m++)
2220  {
2221  Tensordata[R * (Q * (P * (M * i + j) + k) + l) + m] =
2222  source_tensor(i, j, k, l, m);
2223  }
2224  }
2225  }
2226  }
2227  }
2228  }
2229 
2230  /// Copy assignement
2231  RankFiveTensor& operator=(const RankFiveTensor& source_tensor)
2232  {
2233  // Don't create a new matrix if the assignement is the identity
2234  if (this != &source_tensor)
2235  {
2236  // Check row and column length
2237  unsigned long n = source_tensor.nindex1();
2238  unsigned long m = source_tensor.nindex2();
2239  unsigned long p = source_tensor.nindex3();
2240  unsigned long q = source_tensor.nindex4();
2241  unsigned long r = source_tensor.nindex5();
2242  // Resize the tensor to be the same size as the old tensor
2243  if ((N != n) || (M != m) || (P != p) || (Q != q) || (R != r))
2244  {
2245  resize(n, m, p, q, r);
2246  }
2247 
2248  // Copy entries across from the other matrix
2249  for (unsigned long i = 0; i < N; i++)
2250  {
2251  for (unsigned long j = 0; j < M; j++)
2252  {
2253  for (unsigned long k = 0; k < P; k++)
2254  {
2255  for (unsigned long l = 0; l < Q; l++)
2256  {
2257  for (unsigned long m = 0; m < R; m++)
2258  {
2259  (*this)(i, j, k, l, m) = source_tensor(i, j, k, l, m);
2260  }
2261  }
2262  }
2263  }
2264  }
2265  }
2266  // Return reference to object itself (i.e. de-reference this pointer)
2267  return *this;
2268  }
2269 
2270 
2271  /// One parameter constructor produces a nxnxnxnxn tensor
2272  RankFiveTensor(const unsigned long& n)
2273  {
2274  // Set row and column lengths
2275  N = n;
2276  M = n;
2277  P = n;
2278  Q = n;
2279  R = n;
2280  // Assign space for the n rows
2281  Tensordata = new T[N * M * P * Q * R];
2282  // Initialise to zero if required
2283 #ifdef OOMPH_INITIALISE_DENSE_MATRICES
2284  initialise(T(0));
2285 #endif
2286  }
2287 
2288  /// Four parameter constructor, general non-square tensor
2289  RankFiveTensor(const unsigned long& n_index1,
2290  const unsigned long& n_index2,
2291  const unsigned long& n_index3,
2292  const unsigned long& n_index4,
2293  const unsigned long& n_index5)
2294  {
2295  // Set row and column lengths
2296  N = n_index1;
2297  M = n_index2;
2298  P = n_index3;
2299  Q = n_index4;
2300  R = n_index5;
2301  // Assign space for the n rows
2302  Tensordata = new T[N * M * P * Q * R];
2303  // Initialise to zero if required
2304 #ifdef OOMPH_INITIALISE_DENSE_MATRICES
2305  initialise(T(0));
2306 #endif
2307  }
2308 
2309 
2310  /// Four parameter constructor, general non-square tensor
2311  RankFiveTensor(const unsigned long& n_index1,
2312  const unsigned long& n_index2,
2313  const unsigned long& n_index3,
2314  const unsigned long& n_index4,
2315  const unsigned long& n_index5,
2316  const T& initial_val)
2317  {
2318  // Set row and column lengths
2319  N = n_index1;
2320  M = n_index2;
2321  P = n_index3;
2322  Q = n_index4;
2323  R = n_index5;
2324  // Assign space for the n rows
2325  Tensordata = new T[N * M * P * Q * R];
2326  // Initialise to the initial value
2327  initialise(initial_val);
2328  }
2329 
2330  /// Destructor: delete the pointers
2332  {
2333  delete[] Tensordata;
2334  Tensordata = 0;
2335  }
2336 
2337  /// Resize to a square nxnxnxn tensor
2338  void resize(const unsigned long& n)
2339  {
2340  resize(n, n, n, n, n);
2341  }
2342 
2343  /// Resize to a general tensor
2344  void resize(const unsigned long& n_index1,
2345  const unsigned long& n_index2,
2346  const unsigned long& n_index3,
2347  const unsigned long& n_index4,
2348  const unsigned long& n_index5)
2349  {
2350  // If the sizes have not changed do nothing
2351  if ((n_index1 == N) && (n_index2 == M) && (n_index3 == P) &&
2352  (n_index4 == Q) && (n_index5 == R))
2353  {
2354  return;
2355  }
2356  // Store old sizes
2357  unsigned long n_old = N, m_old = M, p_old = P, q_old = Q, r_old = R;
2358  // Reassign the sizes
2359  N = n_index1;
2360  M = n_index2;
2361  P = n_index3;
2362  Q = n_index4;
2363  R = n_index5;
2364  // Store pointer to old matrix data
2365  T* temp_tensor = Tensordata;
2366  // Re-create Tensordata in new size
2367  Tensordata = new T[N * M * P * Q * R];
2368 #ifdef OOMPH_INITIALISE_DENSE_MATRICES
2369  initialise(T(0));
2370 #endif
2371  // Transfer values
2372  unsigned long n_copy, m_copy, p_copy, q_copy, r_copy;
2373  n_copy = std::min(n_old, n_index1);
2374  m_copy = std::min(m_old, n_index2);
2375  p_copy = std::min(p_old, n_index3);
2376  q_copy = std::min(q_old, n_index4);
2377  r_copy = std::min(r_old, n_index5);
2378  // If matrix has values, transfer them to new matrix
2379  // Loop over rows
2380  for (unsigned long i = 0; i < n_copy; i++)
2381  {
2382  // Loop over columns
2383  for (unsigned long j = 0; j < m_copy; j++)
2384  {
2385  // Loop over columns
2386  for (unsigned long k = 0; k < p_copy; k++)
2387  {
2388  // Loop over columns
2389  for (unsigned long l = 0; l < q_copy; l++)
2390  {
2391  // Loop over columns
2392  for (unsigned long m = 0; m < r_copy; m++)
2393  {
2394  // Transfer values from temp_tensor
2395  Tensordata[R * (Q * (M * P * i + P * j + k) + l) + m] =
2396  temp_tensor[r_old *
2397  (q_old * (m_old * p_old * i + p_old * j + k) +
2398  l) +
2399  m];
2400  }
2401  }
2402  }
2403  }
2404  }
2405  // Now kill storage for old tensor
2406  delete[] temp_tensor;
2407  }
2408 
2409  /// Resize to a general tensor
2410  void resize(const unsigned long& n_index1,
2411  const unsigned long& n_index2,
2412  const unsigned long& n_index3,
2413  const unsigned long& n_index4,
2414  const unsigned long& n_index5,
2415  const T& initial_value)
2416  {
2417  // If the sizes have not changed do nothing
2418  if ((n_index1 == N) && (n_index2 == M) && (n_index3 == P) &&
2419  (n_index4 == Q) && (n_index5 == R))
2420  {
2421  return;
2422  }
2423  // Store old sizes
2424  unsigned long n_old = N, m_old = M, p_old = P, q_old = Q, r_old = R;
2425  // Reassign the sizes
2426  N = n_index1;
2427  M = n_index2;
2428  P = n_index3;
2429  Q = n_index4;
2430  R = n_index5;
2431  // Store triple pointer to old matrix data
2432  T* temp_tensor = Tensordata;
2433  // Re-create Tensordata in new size
2434  Tensordata = new T[N * M * P * Q * R];
2435  // Initialise the newly allocated storage
2436  initialise(initial_value);
2437 
2438  // Transfer values
2439  unsigned long n_copy, m_copy, p_copy, q_copy, r_copy;
2440  n_copy = std::min(n_old, n_index1);
2441  m_copy = std::min(m_old, n_index2);
2442  p_copy = std::min(p_old, n_index3);
2443  q_copy = std::min(q_old, n_index4);
2444  r_copy = std::min(r_old, n_index5);
2445  // If matrix has values, transfer them to new matrix
2446  // Loop over rows
2447  for (unsigned long i = 0; i < n_copy; i++)
2448  {
2449  // Loop over columns
2450  for (unsigned long j = 0; j < m_copy; j++)
2451  {
2452  // Loop over columns
2453  for (unsigned long k = 0; k < p_copy; k++)
2454  {
2455  // Loop over columns
2456  for (unsigned long l = 0; l < q_copy; l++)
2457  {
2458  // Loop over columns
2459  for (unsigned long m = 0; m < r_copy; m++)
2460  {
2461  // Transfer values from temp_tensor
2462  Tensordata[R * (Q * (M * P * i + P * j + k) + l) + m] =
2463  temp_tensor[r_old *
2464  (q_old * (m_old * p_old * i + p_old * j + k) +
2465  l) +
2466  m];
2467  }
2468  }
2469  }
2470  }
2471  }
2472  // Now kill storage for old tensor
2473  delete[] temp_tensor;
2474  }
2475 
2476  /// Initialise all values in the tensor to val
2477  void initialise(const T& val)
2478  {
2479  for (unsigned long i = 0; i < (N * M * P * Q * R); ++i)
2480  {
2481  Tensordata[i] = val;
2482  }
2483  }
2484 
2485  /// Return the range of index 1 of the tensor
2486  unsigned long nindex1() const
2487  {
2488  return N;
2489  }
2490 
2491  /// Return the range of index 2 of the tensor
2492  unsigned long nindex2() const
2493  {
2494  return M;
2495  }
2496 
2497  /// Return the range of index 3 of the tensor
2498  unsigned long nindex3() const
2499  {
2500  return P;
2501  }
2502 
2503  /// Return the range of index 4 of the tensor
2504  unsigned long nindex4() const
2505  {
2506  return Q;
2507  }
2508 
2509  /// Return the range of index 5 of the tensor
2510  unsigned long nindex5() const
2511  {
2512  return R;
2513  }
2514 
2515  /// Overload the round brackets to give access as a(i,j,k,l,m)
2516  inline T& operator()(const unsigned long& i,
2517  const unsigned long& j,
2518  const unsigned long& k,
2519  const unsigned long& l,
2520  const unsigned long& m)
2521  {
2522 #ifdef RANGE_CHECKING
2523  this->range_check(i, j, k, l, m);
2524 #endif
2525  return Tensordata[R * (Q * (P * (M * i + j) + k) + l) + m];
2526  }
2527 
2528  /// Overload a const version for read-only access as a(i,j,k,l,m)
2529  inline T operator()(const unsigned long& i,
2530  const unsigned long& j,
2531  const unsigned long& k,
2532  const unsigned long& l,
2533  const unsigned long& m) const
2534  {
2535 #ifdef RANGE_CHECKING
2536  this->range_check(i, j, k, l, m);
2537 #endif
2538  return Tensordata[R * (Q * (P * (M * i + j) + k) + l) + m];
2539  }
2540 
2541  /// Direct access to internal storage of data in flat-packed C-style
2542  /// column-major format. WARNING: Only for experienced users. Only
2543  /// use this if raw speed is of the essence, as in the solid mechanics
2544  /// problems.
2545  inline T& raw_direct_access(const unsigned long& i)
2546  {
2547  return Tensordata[i];
2548  }
2549 
2550 
2551  /// Direct access to internal storage of data in flat-packed C-style
2552  /// column-major format. WARNING: Only for experienced users. Only
2553  /// use this if raw speed is of the essence, as in the solid mechanics
2554  /// problems.
2555  inline const T& raw_direct_access(const unsigned long& i) const
2556  {
2557  return Tensordata[i];
2558  }
2559 
2560  /// Caculate the offset in flat-packed Cy-style, column-major format,
2561  /// required for a given i,j,k. WARNING: Only for experienced users. Only
2562  /// use this if raw speed is of the essence, as in the solid mechanics
2563  /// problems.
2564  unsigned offset(const unsigned long& i,
2565  const unsigned long& j,
2566  const unsigned long& k) const
2567  {
2568  return (R * (Q * (P * (M * i + j) + k) + 0) + 0);
2569  }
2570  };
2571 
2572 
2573  /// ///////////////////////////////////////////////////////////////
2574  /// ///////////////////////////////////////////////////////////////
2575  /// ///////////////////////////////////////////////////////////////
2576 
2577  //=======================================================================
2578  /// A class for compressed column matrices: a sparse matrix format
2579  /// The class is passed as the MATRIX_TYPE paramater so that the base
2580  /// class can use the specific access functions in the round-bracket
2581  /// operator.
2582  //=======================================================================
2583  template<class T>
2584  class CCMatrix : public SparseMatrix<T, CCMatrix<T>>
2585  {
2586  public:
2587  /// Default constructor
2589  {
2590  Row_index = 0;
2591  Column_start = 0;
2592  }
2593 
2594 
2595  /// Constructor: Pass vector of values, vector of row indices,
2596  /// vector of column starts and number of rows (can be suppressed
2597  /// for square matrices). Number of nonzero entries is read
2598  /// off from value, so make sure the vector has been shrunk
2599  /// to its correct length.
2601  const Vector<int>& row_index_,
2602  const Vector<int>& column_start_,
2603  const unsigned long& n,
2604  const unsigned long& m)
2605  : SparseMatrix<T, CCMatrix<T>>()
2606  {
2607  Row_index = 0;
2608  Column_start = 0;
2609  build(value, row_index_, column_start_, n, m);
2610  }
2611 
2612 
2613  /// Copy constructor
2614  CCMatrix(const CCMatrix& source_matrix)
2615  : SparseMatrix<T, CCMatrix<T>>(source_matrix)
2616  {
2617  // NNz, N and M are set the the copy constructor of the SparseMatrix
2618  // called above
2619 
2620  // Row indices stored in C-style array
2621  Row_index = new int[this->Nnz];
2622 
2623  // Assign:
2624  for (unsigned long i = 0; i < this->Nnz; i++)
2625  {
2626  Row_index[i] = source_matrix.row_index()[i];
2627  }
2628 
2629  // Column start:
2630  Column_start = new int[this->M + 1];
2631 
2632  // Assign:
2633  for (unsigned long i = 0; i <= this->M; i++)
2634  {
2635  Column_start[i] = source_matrix.column_start()[i];
2636  }
2637  }
2638 
2639  /// Broken assignment operator
2640  void operator=(const CCMatrix&) = delete;
2641 
2642 
2643  /// Destructor, delete any allocated memory
2644  virtual ~CCMatrix()
2645  {
2646  delete[] Row_index;
2647  Row_index = 0;
2648  delete[] Column_start;
2649  Column_start = 0;
2650  }
2651 
2652  /// Access function that will be called by the read-only
2653  /// round-bracket operator (const)
2654  T get_entry(const unsigned long& i, const unsigned long& j) const
2655  {
2656 #ifdef RANGE_CHECKING
2657  this->range_check(i, j);
2658 #endif
2659  for (long k = Column_start[j]; k < Column_start[j + 1]; k++)
2660  {
2661  if (unsigned(Row_index[k]) == i)
2662  {
2663  return this->Value[k];
2664  }
2665  }
2666  return this->Zero;
2667  }
2668 
2669  /// Read-write access is not permitted for these matrices and is
2670  /// deliberately broken.
2671  T& entry(const unsigned long& i, const unsigned long& j)
2672  {
2673  std::string error_string =
2674  "Non-const access not provided for the CCMatrix<T> class\n";
2675  error_string +=
2676  "It is not possible to use round-bracket access: M(i,j)\n";
2677  error_string += "if M is not declared as const.\n";
2678  error_string += "The solution (albeit ugly) is to create const reference "
2679  "to the matrix\n";
2680  error_string += " const CCMatrix<T>& read_M = M;\n";
2681  error_string += "Then read_M(i,j) is permitted\n";
2682 
2683  throw OomphLibError(
2684  error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2685 
2686  // Dummy return
2687  T dummy;
2688  return dummy;
2689  }
2690 
2691  /// Access to C-style column_start array
2693  {
2694  return Column_start;
2695  }
2696 
2697  /// Access to C-style column_start array (const version)
2698  const int* column_start() const
2699  {
2700  return Column_start;
2701  }
2702 
2703  /// Access to C-style row index array
2704  int* row_index()
2705  {
2706  return Row_index;
2707  }
2708 
2709  /// Access to C-style row index array (const version)
2710  const int* row_index() const
2711  {
2712  return Row_index;
2713  }
2714 
2715  /// Output the "bottom right" entry regardless of it being
2716  /// zero or not (this allows automatic detection of matrix size in
2717  /// e.g. matlab, python).
2718  void output_bottom_right_zero_helper(std::ostream& outfile) const
2719  {
2720  int last_row = this->N - 1;
2721  int last_col_local = this->M - 1;
2722 
2723  // Use this strange thingy because of the CRTP discussed above.
2724  T last_value = this->operator()(last_row, last_col_local);
2725 
2726  if (last_value == T(0))
2727  {
2728  outfile << last_row << " " << last_col_local << " " << T(0)
2729  << std::endl;
2730  }
2731  }
2732 
2733  /// Indexed output function to print a matrix to the
2734  /// stream outfile as i,j,a(i,j) for a(i,j)!=0 only.
2735  void sparse_indexed_output_helper(std::ostream& outfile) const
2736  {
2737  for (unsigned long j = 0; j < this->N; j++)
2738  {
2739  for (long k = Column_start[j]; k < Column_start[j + 1]; k++)
2740  {
2741  outfile << Row_index[k] << " " << j << " " << this->Value[k]
2742  << std::endl;
2743  }
2744  }
2745  }
2746 
2747  /// Wipe matrix data and set all values to 0.
2749 
2750 
2751  /// Build matrix from compressed representation.
2752  /// Number of nonzero entries is read
2753  /// off from value, so make sure the vector has been shrunk
2754  /// to its correct length.
2755  void build(const Vector<T>& value,
2756  const Vector<int>& row_index,
2757  const Vector<int>& column_start,
2758  const unsigned long& n,
2759  const unsigned long& m);
2760 
2761  /// Function to build matrix from pointers to arrays
2762  /// which hold the column starts, row indices and non-zero values.
2763  /// The final parameters specifies the number of rows and columns.
2764  /// Note that, as the name suggests, this function does not
2765  /// make a copy of the data pointed to by the first three arguments!
2767  int* row_index,
2768  int* column_start,
2769  const unsigned long& nnz,
2770  const unsigned long& n,
2771  const unsigned long& m);
2772 
2773 
2774  protected:
2775  /// Row index
2777 
2778  /// Start index for column
2780  };
2781 
2782  /// ////////////////////////////////////////////////////////////////
2783  /// ////////////////////////////////////////////////////////////////
2784  /// ////////////////////////////////////////////////////////////////
2785 
2786 
2787  //=================================================================
2788  /// A class for compressed column matrices that store doubles
2789  //=================================================================
2790  class CCDoubleMatrix : public DoubleMatrixBase, public CCMatrix<double>
2791  {
2792  public:
2793  /// Default constructor
2794  CCDoubleMatrix();
2795 
2796  /// Constructor: Pass vector of values, vector of row indices,
2797  /// vector of column starts and number of rows (can be suppressed
2798  /// for square matrices). Number of nonzero entries is read
2799  /// off from value, so make sure the vector has been shrunk
2800  /// to its correct length.
2802  const Vector<int>& row_index_,
2803  const Vector<int>& column_start_,
2804  const unsigned long& n,
2805  const unsigned long& m);
2806 
2807  /// Broken copy constructor
2808  CCDoubleMatrix(const CCDoubleMatrix& matrix) = delete;
2809 
2810  /// Broken assignment operator
2811  void operator=(const CCDoubleMatrix&) = delete;
2812 
2813  /// Destructor: Kill the LU factors if they have been setup.
2814  virtual ~CCDoubleMatrix();
2815 
2816  /// Return the number of rows of the matrix
2817  inline unsigned long nrow() const
2818  {
2819  return CCMatrix<double>::nrow();
2820  }
2821 
2822  /// Return the number of columns of the matrix
2823  inline unsigned long ncol() const
2824  {
2825  return CCMatrix<double>::ncol();
2826  }
2827 
2828  /// Overload the round-bracket access operator to provide
2829  /// read-only (const) access to the data
2830  inline double operator()(const unsigned long& i,
2831  const unsigned long& j) const
2832  {
2833  return CCMatrix<double>::get_entry(i, j);
2834  }
2835 
2836  /// LU decomposition using SuperLU
2837  virtual void ludecompose();
2838 
2839  /// LU back solve for given RHS
2840  virtual void lubksub(DoubleVector& rhs);
2841 
2842  /// Multiply the matrix by the vector x: soln=Ax
2843  void multiply(const DoubleVector& x, DoubleVector& soln) const;
2844 
2845  /// Multiply the transposed matrix by the vector x: soln=A^T x
2846  void multiply_transpose(const DoubleVector& x, DoubleVector& soln) const;
2847 
2848 
2849  /// Function to multiply this matrix by the CCDoubleMatrix matrix_in
2850  /// The multiplication method used can be selected using the flag
2851  /// Matrix_matrix_multiply_method. By default Method 2 is used.
2852  /// Method 1: First runs through this matrix and matrix_in to find the
2853  /// storage
2854  /// requirements for result - arrays of the correct size are
2855  /// then allocated before performing the calculation.
2856  /// Minimises memory requirements but more costly.
2857  /// Method 2: Grows storage for values and column indices of result 'on the
2858  /// fly' using an array of maps. Faster but more memory
2859  /// intensive.
2860  /// Method 3: Grows storage for values and column indices of result 'on the
2861  /// fly' using a vector of vectors. Not particularly impressive
2862  /// on the platforms we tried...
2863  void multiply(const CCDoubleMatrix& matrix_in, CCDoubleMatrix& result);
2864 
2865 
2866  /// For every row, find the maximum absolute value of the
2867  /// entries in this row. Set all values that are less than alpha times
2868  /// this maximum to zero and return the resulting matrix in
2869  /// reduced_matrix. Note: Diagonal entries are retained regardless
2870  /// of their size.
2871  void matrix_reduction(const double& alpha, CCDoubleMatrix& reduced_matrix);
2872 
2873  /// Access function to Matrix_matrix_multiply_method, the flag
2874  /// which determines the matrix matrix multiplication method used.
2875  /// Method 1: First runs through this matrix and matrix_in to find the
2876  /// storage
2877  /// requirements for result - arrays of the correct size are
2878  /// then allocated before performing the calculation.
2879  /// Minimises memory requirements but more costly.
2880  /// Method 2: Grows storage for values and column indices of result 'on the
2881  /// fly' using an array of maps. Faster but more memory
2882  /// intensive.
2883  /// Method 3: Grows storage for values and column indices of result 'on the
2884  /// fly' using a vector of vectors. Not particularly impressive
2885  /// on the platforms we tried...
2887  {
2889  }
2890 
2891  private:
2892  /// Flag to determine which matrix-matrix multiplication method is used.
2894  };
2895 
2896 
2897  /// //////////////////////////////////////////////////////////////////////
2898  /// //////////////////////////////////////////////////////////////////////
2899  /// //////////////////////////////////////////////////////////////////////
2900 
2901 
2902  //============================================================================
2903  /// Constructor to build a square n by n matrix
2904  //============================================================================
2905  template<class T>
2906  DenseMatrix<T>::DenseMatrix(const unsigned long& n)
2907  {
2908  // Set row and column lengths
2909  N = n;
2910  M = n;
2911  // Assign space for the n rows
2912  Matrixdata = new T[n * n];
2913  // Initialise to zero if required
2914 #ifdef OOMPH_INITIALISE_DENSE_MATRICES
2915  initialise(T(0));
2916 #endif
2917  }
2918 
2919 
2920  //============================================================================
2921  /// Constructor to build a matrix with n rows and m columns
2922  //============================================================================
2923  template<class T>
2924  DenseMatrix<T>::DenseMatrix(const unsigned long& n, const unsigned long& m)
2925  {
2926  // Set row and column lengths
2927  N = n;
2928  M = m;
2929  // Assign space for the n rows
2930  Matrixdata = new T[n * m];
2931 #ifdef OOMPH_INITIALISE_DENSE_MATRICES
2932  initialise(T(0));
2933 #endif
2934  }
2935 
2936  //============================================================================
2937  /// Constructor to build a matrix with n rows and m columns,
2938  /// with initial value initial_val
2939  //============================================================================
2940  template<class T>
2941  DenseMatrix<T>::DenseMatrix(const unsigned long& n,
2942  const unsigned long& m,
2943  const T& initial_val)
2944  {
2945  // Set row and column lengths
2946  N = n;
2947  M = m;
2948  // Assign space for the n rows
2949  Matrixdata = new T[n * m];
2950  initialise(initial_val);
2951  }
2952 
2953 
2954  //============================================================================
2955  /// Resize to a non-square n_row x m_col matrix,
2956  /// where any values already present will be transfered.
2957  //============================================================================
2958  template<class T>
2959  void DenseMatrix<T>::resize(const unsigned long& n, const unsigned long& m)
2960  {
2961  // If the sizes are the same, do nothing
2962  if ((n == N) && (m == M))
2963  {
2964  return;
2965  }
2966  // Store old sizes
2967  unsigned long n_old = N, m_old = M;
2968  // Reassign the sizes
2969  N = n;
2970  M = m;
2971  // Store double pointer to old matrix data
2972  T* temp_matrix = Matrixdata;
2973 
2974  // Re-create Matrixdata in new size
2975  Matrixdata = new T[n * m];
2976  // Initialise to zero
2977 #ifdef OOMPH_INITIALISE_DENSE_MATRICES
2978  initialise(T(0));
2979 #endif
2980 
2981  // Transfer previously existing values
2982  unsigned long n_copy, m_copy;
2983  n_copy = std::min(n_old, n);
2984  m_copy = std::min(m_old, m);
2985 
2986  // If matrix has values, transfer them to new matrix
2987  // Loop over rows
2988  for (unsigned long i = 0; i < n_copy; i++)
2989  {
2990  // Loop over columns
2991  for (unsigned long j = 0; j < m_copy; j++)
2992  {
2993  // Transfer values from temp_matrix
2994  Matrixdata[m * i + j] = temp_matrix[m_old * i + j];
2995  }
2996  }
2997 
2998  // Now kill storage for old matrix
2999  delete[] temp_matrix;
3000  }
3001 
3002 
3003  //============================================================================
3004  /// Resize to a non-square n_row x m_col matrix and initialize the
3005  /// new entries to specified value.
3006  //============================================================================
3007  template<class T>
3008  void DenseMatrix<T>::resize(const unsigned long& n,
3009  const unsigned long& m,
3010  const T& initial_value)
3011  {
3012  // If the size is not changed, just return
3013  if ((n == N) && (m == M))
3014  {
3015  return;
3016  }
3017  // Store old sizes
3018  unsigned long n_old = N, m_old = M;
3019  // Reassign the sizes
3020  N = n;
3021  M = m;
3022  // Store double pointer to old matrix data
3023  T* temp_matrix = Matrixdata;
3024  // Re-create Matrixdata in new size
3025  Matrixdata = new T[n * m];
3026  // Assign initial value (will use the newly allocated data)
3027  initialise(initial_value);
3028 
3029  // Transfering values
3030  unsigned long n_copy, m_copy;
3031  n_copy = std::min(n_old, n);
3032  m_copy = std::min(m_old, m);
3033  // If matrix has values, transfer them to temp_matrix
3034  // Loop over rows
3035  for (unsigned long i = 0; i < n_copy; i++)
3036  {
3037  // Loop over columns
3038  for (unsigned long j = 0; j < m_copy; j++)
3039  {
3040  // Transfer values to temp_matrix
3041  Matrixdata[m * i + j] = temp_matrix[m_old * i + j];
3042  }
3043  }
3044 
3045  // Now kill storage for old matrix
3046  delete[] temp_matrix;
3047  }
3048 
3049 
3050  //============================================================================
3051  /// Output function to print a matrix row-by-row to the stream outfile
3052  //============================================================================
3053  template<class T>
3054  void DenseMatrix<T>::output(std::ostream& outfile) const
3055  {
3056  // Loop over the rows
3057  for (unsigned i = 0; i < N; i++)
3058  {
3059  // Loop over the columne
3060  for (unsigned j = 0; j < M; j++)
3061  {
3062  outfile << (*this)(i, j) << " ";
3063  }
3064  // Put in a newline
3065  outfile << std::endl;
3066  }
3067  }
3068 
3069 
3070  //============================================================================
3071  /// Output function to print a matrix row-by-row to a file. Specify filename.
3072  //============================================================================
3073  template<class T>
3075  {
3076  // Open file
3077  std::ofstream some_file;
3078  some_file.open(filename.c_str());
3079 
3080  output(some_file);
3081  some_file.close();
3082  }
3083 
3084 
3085  //============================================================================
3086  /// Indexed output as i,j,a(i,j)
3087  //============================================================================
3088  template<class T>
3089  void DenseMatrix<T>::indexed_output(std::ostream& outfile) const
3090  {
3091  // Loop over the rows
3092  for (unsigned i = 0; i < N; i++)
3093  {
3094  // Loop over the columns
3095  for (unsigned j = 0; j < M; j++)
3096  {
3097  outfile << i << " " << j << " " << (*this)(i, j) << std::endl;
3098  }
3099  }
3100  }
3101 
3102 
3103  //============================================================================
3104  /// Indexed output function to print a matrix to a
3105  /// file as i,j,a(i,j). Specify filename.
3106  //============================================================================
3107  template<class T>
3109  {
3110  // Open file
3111  std::ofstream some_file;
3112  some_file.open(filename.c_str());
3113  indexed_output(some_file);
3114  some_file.close();
3115  }
3116 
3117 
3118  //============================================================================
3119  /// Output the "bottom right" entry regardless of it being
3120  /// zero or not (this allows automatic detection of matrix size in
3121  /// e.g. matlab, python).
3122  //============================================================================
3123  template<class T>
3125  std::ostream& outfile) const
3126  {
3127  int last_row = this->N - 1;
3128  int last_col = this->M - 1;
3129 
3130  // Use this strange thingy because of the CRTP discussed above.
3131  T last_value = this->operator()(last_row, last_col);
3132 
3133  if (last_value == T(0))
3134  {
3135  outfile << last_row << " " << last_col << " " << T(0) << std::endl;
3136  }
3137  }
3138 
3139  //============================================================================
3140  /// Sparse indexed output as i,j,a(i,j) for a(i,j)!=0 only.
3141  //============================================================================
3142  template<class T>
3143  void DenseMatrix<T>::sparse_indexed_output_helper(std::ostream& outfile) const
3144  {
3145  // Loop over the rows
3146  for (unsigned i = 0; i < N; i++)
3147  {
3148  // Loop over the column
3149  for (unsigned j = 0; j < M; j++)
3150  {
3151  if ((*this)(i, j) != T(0))
3152  {
3153  outfile << i << " " << j << " " << (*this)(i, j) << std::endl;
3154  }
3155  }
3156  }
3157  }
3158 
3159 
3160  /// /////////////////////////////////////////////////////////////////////
3161  /// /////////////////////////////////////////////////////////////////////
3162  /// /////////////////////////////////////////////////////////////////////
3163 
3164 
3165  //=============================================================================
3166  /// Wipe matrix data and set all values to 0.
3167  //=============================================================================
3168  template<class T>
3170  {
3171  // delete any previously allocated storage
3172  if (this->Value != 0)
3173  {
3174  delete[] this->Value;
3175  this->Value = 0;
3176  }
3177  if (this->Row_index != 0)
3178  {
3179  delete[] this->Row_index;
3180  this->Row_index = 0;
3181  }
3182  if (this->Column_start != 0)
3183  {
3184  delete[] this->Column_start;
3185  this->Column_start = 0;
3186  }
3187  this->Nnz = 0;
3188  this->N = 0;
3189  this->M = 0;
3190  }
3191 
3192 
3193  //=============================================================================
3194  /// Build matrix from compressed representation.
3195  /// Note that, as the name suggests, this function does not
3196  /// make a copy of the data pointed to by the first three arguments!
3197  //=============================================================================
3198  template<class T>
3200  int* row_index,
3201  int* column_start,
3202  const unsigned long& nnz,
3203  const unsigned long& n,
3204  const unsigned long& m)
3205  {
3206  // Number of nonzero entries
3207  this->Nnz = nnz;
3208 
3209  // Number of rows
3210  this->N = n;
3211 
3212  // Number of columns
3213  this->M = m;
3214 
3215  // delete any previously allocated storage
3216  if (this->Value != 0)
3217  {
3218  delete[] this->Value;
3219  }
3220  if (this->Row_index != 0)
3221  {
3222  delete[] this->Row_index;
3223  }
3224  if (this->Column_start != 0)
3225  {
3226  delete[] this->Column_start;
3227  }
3228 
3229  // set Value
3230  this->Value = value;
3231 
3232  // set Row_index
3233  this->Row_index = row_index;
3234 
3235  // set Column_start
3236  this->Column_start = column_start;
3237  }
3238 
3239 
3240  //===================================================================
3241  /// Build matrix from compressed representation.
3242  /// Number of nonzero entries is read
3243  /// off from value, so make sure the vector has been shrunk
3244  /// to its correct length.
3245  //===================================================================
3246  template<class T>
3247  void CCMatrix<T>::build(const Vector<T>& value,
3248  const Vector<int>& row_index_,
3249  const Vector<int>& column_start_,
3250  const unsigned long& n,
3251  const unsigned long& m)
3252  {
3253 #ifdef PARANOID
3254  if (value.size() != row_index_.size())
3255  {
3256  std::ostringstream error_message;
3257  error_message << "length of value " << value.size()
3258  << " and row_index vectors " << row_index_.size()
3259  << " should match " << std::endl;
3260 
3261  throw OomphLibError(
3262  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3263  }
3264 #endif
3265 
3266  // Number of nonzero entries
3267  this->Nnz = value.size();
3268 
3269  // Number of rows
3270  this->N = n;
3271 
3272  // Number of columns
3273  this->M = m;
3274 
3275  // We need to delete any previously allocated storage
3276  if (this->Value != 0)
3277  {
3278  delete[] this->Value;
3279  }
3280  if (this->Row_index != 0)
3281  {
3282  delete[] this->Row_index;
3283  }
3284  if (this->Column_start != 0)
3285  {
3286  delete[] this->Column_start;
3287  }
3288 
3289  // Values stored in C-style array
3290  this->Value = new T[this->Nnz];
3291 
3292  // Row indices stored in C-style array
3293  this->Row_index = new int[this->Nnz];
3294 
3295  // Assign:
3296  for (unsigned long i = 0; i < this->Nnz; i++)
3297  {
3298  this->Value[i] = value[i];
3299  this->Row_index[i] = row_index_[i];
3300  }
3301 
3302  // Column start:
3303  // Find the size and aollcate
3304  unsigned long n_column_start = column_start_.size();
3305  this->Column_start = new int[n_column_start];
3306 
3307  // Assign:
3308  for (unsigned long i = 0; i < n_column_start; i++)
3309  {
3310  this->Column_start[i] = column_start_[i];
3311  }
3312  }
3313 
3314  /// //////////////////////////////////////////////////////////////////
3315  /// //////////////////////////////////////////////////////////////////
3316  /// //////////////////////////////////////////////////////////////////
3317 
3318 
3319  //=============================================================================
3320  /// Wipe matrix data and set all values to 0.
3321  //=============================================================================
3322  template<class T>
3324  {
3325  // delete any previously allocated storage
3326  if (this->Value != 0)
3327  {
3328  delete[] this->Value;
3329  this->Value = 0;
3330  }
3331  if (this->Column_index != 0)
3332  {
3333  delete[] this->Column_index;
3334  this->Column_index = 0;
3335  }
3336  if (this->Row_start != 0)
3337  {
3338  delete[] this->Row_start;
3339  this->Row_start = 0;
3340  }
3341  this->Nnz = 0;
3342  this->N = 0;
3343  this->M = 0;
3344  }
3345 
3346 
3347  //=============================================================================
3348  /// Function to build a CRMatrix from pointers to arrays which hold the
3349  /// row starts, column indices and non-zero values
3350  /// Note that, as the name suggests, this function does not
3351  /// make a copy of the data pointed to by the first three arguments!
3352  //=============================================================================
3353  template<class T>
3355  int* column_index_,
3356  int* row_start_,
3357  const unsigned long& nnz,
3358  const unsigned long& n,
3359  const unsigned long& m)
3360  {
3361  // Number of nonzero entries
3362  this->Nnz = nnz;
3363 
3364  // Number of rows
3365  this->N = n;
3366 
3367  // Number of columns
3368  this->M = m;
3369 
3370  // delete any previously allocated storage
3371  if (this->Value != 0)
3372  {
3373  delete[] this->Value;
3374  }
3375  if (this->Column_index != 0)
3376  {
3377  delete[] this->Column_index;
3378  }
3379  if (this->Row_start != 0)
3380  {
3381  delete[] this->Row_start;
3382  }
3383 
3384  // set Value
3385  this->Value = value;
3386 
3387  // set Column_index
3388  this->Column_index = column_index_;
3389 
3390  // set Row_start
3391  this->Row_start = row_start_;
3392  }
3393 
3394 
3395  //=================================================================
3396  /// Build matrix from compressed representation.
3397  /// Number of nonzero entries is read
3398  /// off from value, so make sure the vector has been shrunk
3399  /// to its correct length. The optional final
3400  /// parameter specifies the number of columns. If it is not specified
3401  /// the matrix is assumed to be quadratic.
3402  //=================================================================
3403  template<class T>
3404  void CRMatrix<T>::build(const Vector<T>& value,
3405  const Vector<int>& column_index_,
3406  const Vector<int>& row_start_,
3407  const unsigned long& n,
3408  const unsigned long& m)
3409  {
3410 #ifdef PARANOID
3411  if (value.size() != column_index_.size())
3412  {
3413  std::ostringstream error_message;
3414  error_message << "Must have the same number of values and column indices,"
3415  << "we have " << value.size() << " values and "
3416  << column_index_.size() << " column inidicies."
3417  << std::endl;
3418  throw OomphLibError(
3419  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3420  }
3421 #endif
3422  // Number of nonzero entries
3423  this->Nnz = value.size();
3424 
3425  // Number of rows
3426  this->N = n;
3427 
3428  // Number of columns
3429  this->M = m;
3430 
3431  // We need to delete any previously allocated storage
3432  if (this->Value != 0)
3433  {
3434  delete[] this->Value;
3435  }
3436  if (this->Column_index != 0)
3437  {
3438  delete[] this->Column_index;
3439  }
3440  if (this->Row_start != 0)
3441  {
3442  delete[] this->Row_start;
3443  }
3444 
3445  // Values stored in C-style array
3446  this->Value = new T[this->Nnz];
3447 
3448  // Column indices stored in C-style array
3449  this->Column_index = new int[this->Nnz];
3450 
3451  // Assign:
3452  for (unsigned long i = 0; i < this->Nnz; i++)
3453  {
3454  this->Value[i] = value[i];
3455  this->Column_index[i] = column_index_[i];
3456  }
3457 
3458  // Row start:
3459  // Find the size and allocate
3460  unsigned long n_row_start = row_start_.size();
3461  this->Row_start = new int[n_row_start];
3462 
3463  // Assign:
3464  for (unsigned long i = 0; i < n_row_start; i++)
3465  {
3466  this->Row_start[i] = row_start_[i];
3467  }
3468  }
3469 
3470 
3471  //=================================================================
3472  /// Dummy zero
3473  //=================================================================
3474  template<class T, class MATRIX_TYPE>
3476 
3477 
3478  namespace RRR
3479  {
3480  extern std::string RayStr;
3481  extern bool RayBool;
3482  } // namespace RRR
3483 
3484  //=================================================================
3485  /// Namespace for helper functions for CRDoubleMatrices
3486  //=================================================================
3487  namespace CRDoubleMatrixHelpers
3488  {
3489  /// Create a deep copy of the matrix pointed to by in_matrix_pt
3490  inline void deep_copy(const CRDoubleMatrix* const in_matrix_pt,
3491  CRDoubleMatrix& out_matrix)
3492  {
3493 #ifdef PARANOID
3494  // Is the out matrix built? We need an empty out matrix!
3495  if (out_matrix.built())
3496  {
3497  std::ostringstream err_msg;
3498  err_msg << "The result matrix has been built.\n"
3499  << "Please clear the matrix.\n";
3500  throw OomphLibError(
3501  err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3502  }
3503 
3504  // Check that the in matrix pointer is not null.
3505  if (in_matrix_pt == 0)
3506  {
3507  std::ostringstream err_msg;
3508  err_msg << "The in_matrix_pt is null.\n";
3509  throw OomphLibError(
3510  err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3511  }
3512 
3513  // Check that the in matrix is built.
3514  if (!in_matrix_pt->built())
3515  {
3516  std::ostringstream err_msg;
3517  err_msg << "The in_matrix_pt is null.\n";
3518  throw OomphLibError(
3519  err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3520  }
3521 #endif
3522 
3523  // First set the matrix matrix multiply methods (for both serial and
3524  // distributed)
3526  in_matrix_pt->serial_matrix_matrix_multiply_method();
3527 
3530 
3531 
3532  // The local nrow and nnz of the in matrix
3533  const unsigned in_nrow_local = in_matrix_pt->nrow_local();
3534  const unsigned long in_nnz = in_matrix_pt->nnz();
3535 
3536  // Storage for the values, column indices and row start
3537  double* out_values = new double[in_nnz];
3538  int* out_column_indices = new int[in_nnz];
3539  int* out_row_start = new int[in_nrow_local + 1];
3540 
3541  // The data to copy over
3542  const double* const in_values = in_matrix_pt->value();
3543  const int* const in_column_indices = in_matrix_pt->column_index();
3544  const int* const in_row_start = in_matrix_pt->row_start();
3545 
3546  // Copy the data
3547  std::copy(in_values, in_values + in_nnz, out_values);
3548 
3549  std::copy(
3550  in_column_indices, in_column_indices + in_nnz, out_column_indices);
3551 
3552  std::copy(
3553  in_row_start, in_row_start + (in_nrow_local + 1), out_row_start);
3554 
3555  // Build the matrix
3556  out_matrix.build(in_matrix_pt->distribution_pt());
3557 
3558  out_matrix.build_without_copy(in_matrix_pt->ncol(),
3559  in_nnz,
3560  out_values,
3561  out_column_indices,
3562  out_row_start);
3563 
3564  // The only thing we haven't copied over is the default linear solver
3565  // pointer, but I cannot figure out how to copy over a solver since
3566  // I do not know what it is.
3567  } // EoFunc deep_copy
3568 
3569  /// Builds a uniformly distributed matrix.
3570  /// A locally replicated matrix is constructed then redistributed using
3571  /// OOMPH-LIB's default uniform row distribution.
3572  /// This is memory intensive thus should be used for
3573  /// testing or small problems only.
3574  /// The resulting matrix (mat_out) must not have been built.
3576  const unsigned& nrow,
3577  const unsigned& ncol,
3578  const OomphCommunicator* const comm_pt,
3579  const Vector<double>& values,
3580  const Vector<int>& column_indicies,
3581  const Vector<int>& row_start,
3582  CRDoubleMatrix& mat_out);
3583 
3584 
3585  /// Calculates the infinity (maximum) norm of a DenseMartrix of
3586  /// CRDoubleMatrices as if it was one large matrix.
3587  /// This avoids creating a concatenation of the sub-blocks just to calculate
3588  /// the infinity norm.
3589  double inf_norm(const DenseMatrix<CRDoubleMatrix*>& matrix_pt);
3590 
3591  /// Calculates the largest Gershgorin disc whilst preserving the sign. Let
3592  /// A be an n by n matrix, with entries aij. For \f$ i \in \{ 1,...,n \} \f$
3593  /// let \f$ R_i = \sum_{i\neq j} |a_{ij}| \f$ be the sum of the absolute
3594  /// values of the non-diagonal entries in the i-th row. Let \f$ D(a_{ii},R_i) \f$
3595  /// be the closed disc centered at \f$ a_{ii} \f$ with radius \f$ R_i \f$,
3596  /// such a disc is called a Gershgorin disc.
3597  ///
3598  /// \n
3599  ///
3600  /// We calculate \f$ |D(a_{ii},R_i)|_{max} \f$and multiply by the sign of
3601  /// the diagonal entry.
3602  ///
3603  /// \n
3604  ///
3605  /// The DenseMatrix of CRDoubleMatrices are treated as if they are one
3606  /// large matrix. Therefore the dimensions of the sub matrices has to
3607  /// "make sense", there is a paranoid check for this.
3609  const DenseMatrix<CRDoubleMatrix*>& matrix_pt);
3610 
3611  /// Concatenate CRDoubleMatrix matrices.
3612  /// The in matrices are concatenated such that the block structure of the
3613  /// in matrices are preserved in the result matrix. Communication between
3614  /// processors is required. If the block structure of the sub matrices does
3615  /// not need to be preserved, consider using
3616  /// CRDoubleMatrixHelpers::concatenate_without_communication(...).
3617  ///
3618  /// The matrix manipulation functions
3619  /// CRDoubleMatrixHelpers::concatenate(...) and
3620  /// CRDoubleMatrixHelpers::concatenate_without_communication(...)
3621  /// are analogous to the Vector manipulation functions
3622  /// DoubleVectorHelpers::concatenate(...) and
3623  /// DoubleVectorHelpers::concatenate_without_communication(...).
3624  /// Please look at the DoubleVector functions for an illustration of the
3625  /// differences between concatenate(...) and
3626  /// concatenate_without_communication(...).
3627  ///
3628  /// Distribution of the result matrix:
3629  /// If the result matrix does not have a distribution built, then it will be
3630  /// given a uniform row distribution. Otherwise we use the existing
3631  /// distribution. This gives the user the ability to define their own
3632  /// distribution, or save computing power if a distribution has
3633  /// been pre-built.
3634  ///
3635  /// NOTE: ALL the matrices pointed to by matrix_pt has to be built. This is
3636  /// not the case with concatenate_without_communication(...)
3637  void concatenate(const DenseMatrix<CRDoubleMatrix*>& matrix_pt,
3638  CRDoubleMatrix& result_matrix);
3639 
3640  /// Concatenate CRDoubleMatrix matrices.
3641  ///
3642  /// The Vector row_distribution_pt contains the LinearAlgebraDistribution
3643  /// of each block row.
3644  /// The Vector col_distribution_pt contains the LinearAlgebraDistribution
3645  /// of each block column.
3646  /// The DenseMatrix matrix_pt contains pointers to the CRDoubleMatrices
3647  /// to concatenate.
3648  /// The CRDoubleMatrix result_matrix is the result matrix.
3649  ///
3650  /// The result matrix is a permutation of the sub matrices such that the
3651  /// data stays on the same processor when the result matrix is built, there
3652  /// is no communication between processors. Thus the block structure of the
3653  /// sub matrices are NOT preserved in the result matrix. The rows are
3654  /// block-permuted, defined by the concatenation of the distributions in
3655  /// row_distribution_pt. Similarly, the columns are block-permuted, defined
3656  /// by the concatenation of the distributions in col_distribution_pt. For
3657  /// more details on the block-permutation, see
3658  /// LinearAlgebraDistributionHelpers::concatenate(...).
3659  ///
3660  /// If one wishes to preserve the block structure of the sub matrices in the
3661  /// result matrix, consider using CRDoubleMatrixHelpers::concatenate(...),
3662  /// which uses communication between processors to ensure that the block
3663  /// structure of the sub matrices are preserved.
3664  ///
3665  /// The matrix manipulation functions
3666  /// CRDoubleMatrixHelpers::concatenate(...) and
3667  /// CRDoubleMatrixHelpers::concatenate_without_communication(...)
3668  /// are analogous to the Vector manipulation functions
3669  /// DoubleVectorHelpers::concatenate(...) and
3670  /// DoubleVectorHelpers::concatenate_without_communication(...).
3671  /// Please look at the DoubleVector functions for an illustration of the
3672  /// differences between concatenate(...) and
3673  /// concatenate_without_communication(...).
3674  ///
3675  /// Distribution of the result matrix:
3676  /// If the result matrix does not have a distribution built, then it will be
3677  /// given a distribution built from the concatenation of the distributions
3678  /// from row_distribution_pt, see
3679  /// LinearAlgebraDistributionHelpers::concatenate(...) for more detail.
3680  /// Otherwise we use the existing distribution.
3681  /// If there is an existing distribution then it must be the same as the
3682  /// distribution from the concatenation of row distributions as described
3683  /// above.
3684  /// Why don't we always compute the distribution "on the fly"?
3685  /// Because a non-uniform distribution requires communication.
3686  /// All block preconditioner distributions are concatenations of the
3687  /// distributions of the individual blocks.
3689  const Vector<LinearAlgebraDistribution*>& row_distribution_pt,
3690  const Vector<LinearAlgebraDistribution*>& col_distribution_pt,
3691  const DenseMatrix<CRDoubleMatrix*>& matrix_pt,
3692  CRDoubleMatrix& result_matrix);
3693 
3694  /// Concatenate CRDoubleMatrix matrices.
3695  /// This calls the other concatenate_without_communication(...) function,
3696  /// passing block_distribution_pt as both the row_distribution_pt and
3697  /// col_distribution_pt. This should only be called for block square
3698  /// matrices.
3700  const Vector<LinearAlgebraDistribution*>& block_distribution_pt,
3701  const DenseMatrix<CRDoubleMatrix*>& matrix_pt,
3702  CRDoubleMatrix& result_matrix);
3703 
3704  } // namespace CRDoubleMatrixHelpers
3705 
3706 } // namespace oomph
3707 #endif
cstr elem_len * i
Definition: cfortran.h:603
//////////////////////////////////////////////////////////////// ////////////////////////////////////...
Definition: matrices.h:2791
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.
Definition: matrices.h:2830
void operator=(const CCDoubleMatrix &)=delete
Broken assignment operator.
void multiply_transpose(const DoubleVector &x, DoubleVector &soln) const
Multiply the transposed matrix by the vector x: soln=A^T x.
Definition: matrices.cc:715
virtual void lubksub(DoubleVector &rhs)
LU back solve for given RHS.
Definition: matrices.cc:614
void matrix_reduction(const double &alpha, CCDoubleMatrix &reduced_matrix)
For every row, find the maximum absolute value of the entries in this row. Set all values that are le...
Definition: matrices.cc:1149
virtual ~CCDoubleMatrix()
Destructor: Kill the LU factors if they have been setup.
Definition: matrices.cc:597
CCDoubleMatrix(const CCDoubleMatrix &matrix)=delete
Broken copy constructor.
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:2823
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:2817
CCDoubleMatrix()
Default constructor.
Definition: matrices.cc:572
unsigned & matrix_matrix_multiply_method()
Access function to Matrix_matrix_multiply_method, the flag which determines the matrix matrix multipl...
Definition: matrices.h:2886
virtual void ludecompose()
LU decomposition using SuperLU.
Definition: matrices.cc:606
void multiply(const DoubleVector &x, DoubleVector &soln) const
Multiply the matrix by the vector x: soln=Ax.
Definition: matrices.cc:622
unsigned Matrix_matrix_multiply_method
Flag to determine which matrix-matrix multiplication method is used.
Definition: matrices.h:2893
/////////////////////////////////////////////////////////////// /////////////////////////////////////...
Definition: matrices.h:2585
CCMatrix()
Default constructor.
Definition: matrices.h:2588
void sparse_indexed_output_helper(std::ostream &outfile) const
Indexed output function to print a matrix to the stream outfile as i,j,a(i,j) for a(i,...
Definition: matrices.h:2735
CCMatrix(const Vector< T > &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...
Definition: matrices.h:2600
void build_without_copy(T *value, int *row_index, int *column_start, const unsigned long &nnz, const unsigned long &n, const unsigned long &m)
Function to build matrix from pointers to arrays which hold the column starts, row indices and non-ze...
Definition: matrices.h:3199
virtual ~CCMatrix()
Destructor, delete any allocated memory.
Definition: matrices.h:2644
int * Column_start
Start index for column.
Definition: matrices.h:2779
int * Row_index
Row index.
Definition: matrices.h:2776
CCMatrix(const CCMatrix &source_matrix)
Copy constructor.
Definition: matrices.h:2614
T & entry(const unsigned long &i, const unsigned long &j)
Read-write access is not permitted for these matrices and is deliberately broken.
Definition: matrices.h:2671
int * column_start()
Access to C-style column_start array.
Definition: matrices.h:2692
const int * column_start() const
Access to C-style column_start array (const version)
Definition: matrices.h:2698
T 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
void build(const Vector< T > &value, const Vector< int > &row_index, const Vector< int > &column_start, const unsigned long &n, const unsigned long &m)
Build matrix from compressed representation. Number of nonzero entries is read off from value,...
Definition: matrices.h:3247
void operator=(const CCMatrix &)=delete
Broken assignment operator.
const int * row_index() const
Access to C-style row index array (const version)
Definition: matrices.h:2710
void output_bottom_right_zero_helper(std::ostream &outfile) const
Output the "bottom right" entry regardless of it being zero or not (this allows automatic detection o...
Definition: matrices.h:2718
int * row_index()
Access to C-style row index array.
Definition: matrices.h:2704
void clean_up_memory()
Wipe matrix data and set all values to 0.
Definition: matrices.h:3169
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
void sort_entries()
Sorts the entries associated with each row of the matrix in the column index vector and the value vec...
Definition: matrices.cc:1449
int * column_index()
Access to C-style column index array.
Definition: matrices.h:1072
int * row_start()
Access to C-style row_start array.
Definition: matrices.h:1060
virtual ~CRDoubleMatrix()
Destructor.
Definition: matrices.cc:1343
const double * value() const
Access to C-style value array (const version)
Definition: matrices.h:1090
void sparse_indexed_output_with_offset(std::string filename)
Indexed output function to print a matrix to a file as i,j,a(i,j) for a(i,j)!=0 only....
Definition: matrices.h:1031
unsigned & distributed_matrix_matrix_multiply_method()
Access function to Distributed_matrix_matrix_multiply_method, the flag which determines the matrix ma...
Definition: matrices.h:1191
struct oomph::CRDoubleMatrix::CRDoubleMatrixComparisonHelper Comparison_struct
void matrix_reduction(const double &alpha, CRDoubleMatrix &reduced_matrix)
For every row, find the maximum absolute value of the entries in this row. Set all values that are le...
Definition: matrices.cc:2365
void operator=(const CRDoubleMatrix &)=delete
Broken assignment operator.
void multiply_transpose(const DoubleVector &x, DoubleVector &soln) const
Multiply the transposed matrix by the vector x: soln=A^T x.
Definition: matrices.cc:1882
void sparse_indexed_output_helper(std::ostream &outfile) const
Indexed output function to print a matrix to the stream outfile as i,j,a(i,j) for a(i,...
Definition: matrices.h:1023
unsigned Distributed_matrix_matrix_multiply_method
Flag to determine which matrix-matrix multiplication method is used (for distributed matrices)
Definition: matrices.h:1246
virtual void ludecompose()
LU decomposition using SuperLU if matrix is not distributed or distributed onto a single processor.
Definition: matrices.cc:1728
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:1008
void multiply(const DoubleVector &x, DoubleVector &soln) const
Multiply the matrix by the vector x: soln=Ax.
Definition: matrices.cc:1782
void add(const CRDoubleMatrix &matrix_in, CRDoubleMatrix &result_matrix) const
element-wise addition of this matrix with matrix_in.
Definition: matrices.cc:3515
const unsigned & distributed_matrix_matrix_multiply_method() const
Read only access function (const version) to Distributed_matrix_matrix_multiply_method,...
Definition: matrices.h:1202
bool Built
Flag to indicate whether the matrix has been built - i.e. the distribution has been setup AND the mat...
Definition: matrices.h:1253
Vector< int > Index_of_diagonal_entries
Vector whose i'th entry contains the index of the last entry below or on the diagonal of the i'th row...
Definition: matrices.h:1238
double inf_norm() const
returns the inf-norm of this matrix
Definition: matrices.cc:3412
const Vector< int > get_index_of_diagonal_entries() const
Access function: returns the vector Index_of_diagonal_entries. The i-th entry of the vector contains ...
Definition: matrices.h:920
void get_matrix_transpose(CRDoubleMatrix *result) const
Returns the transpose of this matrix.
Definition: matrices.cc:3271
const int * column_index() const
Access to C-style column index array (const version)
Definition: matrices.h:1078
unsigned long nnz() const
Return the number of nonzero entries (the local nnz)
Definition: matrices.h:1096
CRDoubleMatrix * global_matrix() const
if this matrix is distributed then a the equivalent global matrix is built using new and returned....
Definition: matrices.cc:2431
void redistribute(const LinearAlgebraDistribution *const &dist_pt)
The contents of the matrix are redistributed to match the new distribution. In a non-MPI build this m...
Definition: matrices.cc:2575
const int * row_start() const
Access to C-style row_start array (const version)
Definition: matrices.h:1066
bool built() const
access function to the Built flag - indicates whether the matrix has been build - i....
Definition: matrices.h:1210
unsigned & serial_matrix_matrix_multiply_method()
Access function to Serial_matrix_matrix_multiply_method, the flag which determines the matrix matrix ...
Definition: matrices.h:1159
void build_without_copy(const unsigned &ncol, const unsigned &nnz, double *value, int *column_index, int *row_start)
keeps the existing distribution and just matrix that is stored without copying the matrix data
Definition: matrices.cc:1710
void output_bottom_right_zero_helper(std::ostream &outfile) const
Output the "bottom right" entry regardless of it being zero or not (this allows automatic detection o...
Definition: matrices.h:1016
double * value()
Access to C-style value array.
Definition: matrices.h:1084
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:1002
unsigned Serial_matrix_matrix_multiply_method
Flag to determine which matrix-matrix multiplication method is used (for serial (or global) matrices)
Definition: matrices.h:1242
Vector< double > diagonal_entries() const
returns a Vector of diagonal entries of this matrix. This only works with square matrices....
Definition: matrices.cc:3465
CRDoubleMatrix()
Default constructor.
Definition: matrices.cc:1214
bool entries_are_sorted(const bool &doc_unordered_entries=false) const
Runs through the column index vector and checks if the entries follow the regular lexicographical ord...
Definition: matrices.cc:1366
const unsigned & serial_matrix_matrix_multiply_method() const
Read only access function (const version) to Serial_matrix_matrix_multiply_method,...
Definition: matrices.h:1181
CRMatrix< double > CR_matrix
Storage for the Matrix in CR Format.
Definition: matrices.h:1249
virtual void lubksub(DoubleVector &rhs)
LU back solve for given RHS.
Definition: matrices.cc:1749
void build(const LinearAlgebraDistribution *distribution_pt, const unsigned &ncol, const Vector< double > &value, const Vector< int > &column_index, const Vector< int > &row_start)
build method: vector of values, vector of column indices, vector of row starts and number of rows and...
Definition: matrices.cc:1672
double operator()(const unsigned long &i, const unsigned long &j) const
Overload the round-bracket access operator for read-only access. In a distributed matrix i refers to ...
Definition: matrices.h:1053
A class for compressed row matrices, a sparse storage format Once again the recursive template trick ...
Definition: matrices.h:682
CRMatrix()
Default constructor.
Definition: matrices.h:685
void clean_up_memory()
Wipe matrix data and set all values to 0.
Definition: matrices.h:3323
const int * column_index() const
Access to C-style column index array (const version)
Definition: matrices.h:803
CRMatrix(const CRMatrix &source_matrix)
Copy constructor.
Definition: matrices.h:710
const int * row_start() const
Access to C-style row_start array (const version)
Definition: matrices.h:791
T & entry(const unsigned long &i, const unsigned long &j)
The read-write access function is deliberately broken.
Definition: matrices.h:764
int * row_start()
Access to C-style row_start array.
Definition: matrices.h:785
int * Row_start
Start index for row.
Definition: matrices.h:873
int * Column_index
Column index.
Definition: matrices.h:870
void build(const Vector< T > &value, const Vector< int > &column_index, const Vector< int > &row_start, const unsigned long &n, const unsigned long &m)
Build matrix from compressed representation. Number of nonzero entries is read off from value,...
Definition: matrices.h:3404
void operator=(const CRMatrix &)=delete
Broken assignment operator.
int * column_index()
Access to C-style column index array.
Definition: matrices.h:797
void output_bottom_right_zero_helper(std::ostream &outfile) const
Output the "bottom right" entry regardless of it being zero or not (this allows automatic detection o...
Definition: matrices.h:811
virtual ~CRMatrix()
Destructor, delete any allocated memory.
Definition: matrices.h:738
void build_without_copy(T *value, int *column_index, int *row_start, const unsigned long &nnz, const unsigned long &n, const unsigned long &m)
Function to build matrix from pointers to arrays which hold the row starts, column indices and non-ze...
Definition: matrices.h:3354
void sparse_indexed_output_helper(std::ostream &outfile) const
Indexed output function to print a matrix to the stream outfile as i,j,a(i,j) for a(i,...
Definition: matrices.h:828
T 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
CRMatrix(const Vector< T > &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 rows...
Definition: matrices.h:697
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1271
DenseDoubleMatrix()
Constructor, set the default linear solver.
Definition: matrices.cc:139
DenseDoubleMatrix(const DenseDoubleMatrix &matrix)=delete
Broken copy constructor.
virtual void lubksub(DoubleVector &rhs)
LU backsubstitution.
Definition: matrices.cc:202
virtual ~DenseDoubleMatrix()
Destructor.
Definition: matrices.cc:182
void matrix_reduction(const double &alpha, DenseDoubleMatrix &reduced_matrix)
For every row, find the maximum absolute value of the entries in this row. Set all values that are le...
Definition: matrices.cc:481
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.
Definition: matrices.h:1316
virtual void ludecompose()
LU decomposition using DenseLU (default linea solver)
Definition: matrices.cc:192
void multiply_transpose(const DoubleVector &x, DoubleVector &soln) const
Multiply the transposed matrix by the vector x: soln=A^T x.
Definition: matrices.cc:385
void multiply(const DoubleVector &x, DoubleVector &soln) const
Multiply the matrix by the vector x: soln=Ax.
Definition: matrices.cc:294
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:1295
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.
Definition: matrices.h:1308
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:1301
void eigenvalues_by_jacobi(Vector< double > &eigen_val, DenseMatrix< double > &eigen_vect) const
Determine eigenvalues and eigenvectors, using Jacobi rotations. Only for symmetric matrices....
Definition: matrices.cc:224
void operator=(const DenseDoubleMatrix &)=delete
Broken assignment operator.
Dense LU decomposition-based solve of full assembled linear system. VERY inefficient but useful to il...
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
Definition: matrices.h:386
DenseMatrix & operator=(const DenseMatrix &source_matrix)
Copy assignment.
Definition: matrices.h:420
T & 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
void output(std::ostream &outfile) const
Output function to print a matrix row-by-row to the stream outfile.
Definition: matrices.h:3054
T 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
void resize(const unsigned long &n, const unsigned long &m)
Resize to a non-square n x m matrix; any values already present will be transfered.
Definition: matrices.h:2959
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:485
DenseMatrix(const DenseMatrix &source_matrix)
Copy constructor: Deep copy!
Definition: matrices.h:402
DenseMatrix(const unsigned long &n)
Constructor to build a square n by n matrix.
Definition: matrices.h:2906
unsigned long N
Number of rows.
Definition: matrices.h:392
void indexed_output(std::ostream &outfile) const
Indexed output function to print a matrix to the stream outfile as i,j,a(i,j)
Definition: matrices.h:3089
void output(std::string filename) const
Output function to print a matrix row-by-row to a file. Specify filename.
Definition: matrices.h:3074
void indexed_output(std::string filename) const
Indexed output function to print a matrix to a file as i,j,a(i,j). Specify filename.
Definition: matrices.h:3108
DenseMatrix(const unsigned long &n, const unsigned long &m)
Constructor to build a matrix with n rows and m columns.
Definition: matrices.h:2924
virtual ~DenseMatrix()
Destructor, clean up the matrix data.
Definition: matrices.h:478
void sparse_indexed_output_helper(std::ostream &outfile) const
Indexed output function to print a matrix to the stream outfile as i,j,a(i,j) for a(i,...
Definition: matrices.h:3143
DenseMatrix(const unsigned long &n, const unsigned long &m, const T &initial_val)
Constructor to build a matrix with n rows and m columns, with initial value initial_val.
Definition: matrices.h:2941
T * Matrixdata
Internal representation of matrix as a pointer to data.
Definition: matrices.h:389
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:491
void output_bottom_right_zero_helper(std::ostream &outfile) const
Output the "bottom right" entry regardless of it being zero or not (this allows automatic detection o...
Definition: matrices.h:3124
unsigned long M
Number of columns.
Definition: matrices.h:395
void initialise(const T &val)
Initialize all values in the matrix to val.
Definition: matrices.h:514
void resize(const unsigned long &n)
Resize to a square nxn matrix; any values already present will be transfered.
Definition: matrices.h:498
void resize(const unsigned long &n, const unsigned long &m, const T &initial_value)
Resize to a non-square n x m matrix and initialize the new values to initial_value.
Definition: matrices.h:3008
DenseMatrix()
Empty constructor, simply assign the lengths N and M to 0.
Definition: matrices.h:399
Base class for any linear algebra object that is distributable. Just contains storage for the LinearA...
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
unsigned nrow() const
access function to the number of global rows.
unsigned nrow_local() const
access function for the num of local rows on this processor.
unsigned first_row() const
access function for the first row on this processor
Abstract base class for matrices of doubles – adds abstract interfaces for solving,...
Definition: matrices.h:261
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
Definition: matrices.h:296
virtual 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 unsigned long ncol() const =0
Return the number of columns of the matrix.
DoubleMatrixBase(const DoubleMatrixBase &matrix)=delete
Broken copy constructor.
void solve(DoubleVector &rhs)
Complete LU solve (replaces matrix by its LU decomposition and overwrites RHS with solution)....
Definition: matrices.cc:50
virtual double max_residual(const DoubleVector &x, const DoubleVector &rhs)
Find the maximum residual r=b-Ax – generic version, can be overloaded for specific derived classes wh...
Definition: matrices.h:348
LinearSolver * Default_linear_solver_pt
Definition: matrices.h:267
virtual void multiply(const DoubleVector &x, DoubleVector &soln) const =0
Multiply the matrix by the vector x: soln=Ax.
virtual ~DoubleMatrixBase()
virtual (empty) destructor
Definition: matrices.h:286
virtual void multiply_transpose(const DoubleVector &x, DoubleVector &soln) const =0
Multiply the transposed matrix by the vector x: soln=A^T x.
LinearSolver * Linear_solver_pt
Definition: matrices.h:264
virtual void residual(const DoubleVector &x, const DoubleVector &b, DoubleVector &residual_)
Find the residual, i.e. r=b-Ax the residual.
Definition: matrices.h:326
DoubleMatrixBase()
(Empty) constructor.
Definition: matrices.h:271
void operator=(const DoubleMatrixBase &)=delete
Broken assignment operator.
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
LinearSolver *const & linear_solver_pt() const
Return a pointer to the linear solver object (const version)
Definition: matrices.h:302
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
double max() const
returns the maximum coefficient
double * values_pt()
access function to the underlying values
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
unsigned first_row() const
access function for the first row on this processor. If not distributed then this is just zero.
Base class for all linear solvers. This merely defines standard interfaces for linear solvers,...
Definition: linear_solver.h:68
Abstract base class for matrices, templated by the type of object that is stored in them and the type...
Definition: matrices.h:74
virtual void output_bottom_right_zero_helper(std::ostream &outfile) const =0
Output the "bottom right" entry regardless of it being zero or not (this allows automatic detection o...
T & operator()(const unsigned long &i, const unsigned long &j)
Round brackets to give access as a(i,j) for read-write access. The function uses the MATRIX_TYPE temp...
Definition: matrices.h:140
T operator()(const unsigned long &i, const unsigned long &j) const
Round brackets to give access as a(i,j) for read only (we're not providing a general interface for co...
Definition: matrices.h:128
void range_check(const unsigned long &i, const unsigned long &j) const
Range check to catch when an index is out of bounds, if so, it issues a warning message and dies by t...
Definition: matrices.h:78
virtual ~Matrix()
Virtual (empty) destructor.
Definition: matrices.h:114
void sparse_indexed_output(std::ostream &outfile, const unsigned &precision=0, const bool &output_bottom_right_zero=false) const
Indexed output function to print a matrix to the stream outfile as i,j,a(i,j) for a(i,...
Definition: matrices.h:182
virtual void sparse_indexed_output_helper(std::ostream &outfile) const =0
Indexed output function to print a matrix to the stream outfile as i,j,a(i,j) for a(i,...
void operator=(const Matrix &)=delete
Broken assignment operator.
Matrix()
(Empty) constructor
Definition: matrices.h:105
virtual unsigned long nrow() const =0
Return the number of rows of the matrix.
Matrix(const Matrix &matrix)=delete
Broken copy constructor.
void sparse_indexed_output(std::string filename, const unsigned &precision=0, const bool &output_bottom_right_zero=false) const
Indexed output function to print a matrix to the file named filename as i,j,a(i,j) for a(i,...
Definition: matrices.h:226
virtual void output(std::ostream &outfile) const
Output function to print a matrix row-by-row, in the form a(0,0) a(0,1) ... a(1,0) a(1,...
Definition: matrices.h:152
virtual unsigned long ncol() const =0
Return the number of columns of the matrix.
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:54
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
/////////////////////////////////////////////////////////////// /////////////////////////////////////...
Definition: matrices.h:2113
unsigned R
5th Tensor dimension
Definition: matrices.h:2131
T operator()(const unsigned long &i, const unsigned long &j, const unsigned long &k, const unsigned long &l, const unsigned long &m) const
Overload a const version for read-only access as a(i,j,k,l,m)
Definition: matrices.h:2529
const T & raw_direct_access(const unsigned long &i) const
Direct access to internal storage of data in flat-packed C-style column-major format....
Definition: matrices.h:2555
RankFiveTensor(const RankFiveTensor &source_tensor)
Copy constructor: Deep copy.
Definition: matrices.h:2198
unsigned M
2nd Tensor dimension
Definition: matrices.h:2122
unsigned long nindex4() const
Return the range of index 4 of the tensor.
Definition: matrices.h:2504
RankFiveTensor(const unsigned long &n)
One parameter constructor produces a nxnxnxnxn tensor.
Definition: matrices.h:2272
RankFiveTensor()
Empty constructor.
Definition: matrices.h:2195
RankFiveTensor & operator=(const RankFiveTensor &source_tensor)
Copy assignement.
Definition: matrices.h:2231
T & operator()(const unsigned long &i, const unsigned long &j, const unsigned long &k, const unsigned long &l, const unsigned long &m)
Overload the round brackets to give access as a(i,j,k,l,m)
Definition: matrices.h:2516
unsigned Q
4th Tensor dimension
Definition: matrices.h:2128
unsigned long nindex3() const
Return the range of index 3 of the tensor.
Definition: matrices.h:2498
RankFiveTensor(const unsigned long &n_index1, const unsigned long &n_index2, const unsigned long &n_index3, const unsigned long &n_index4, const unsigned long &n_index5, const T &initial_val)
Four parameter constructor, general non-square tensor.
Definition: matrices.h:2311
void resize(const unsigned long &n_index1, const unsigned long &n_index2, const unsigned long &n_index3, const unsigned long &n_index4, const unsigned long &n_index5, const T &initial_value)
Resize to a general tensor.
Definition: matrices.h:2410
void resize(const unsigned long &n)
Resize to a square nxnxnxn tensor.
Definition: matrices.h:2338
void range_check(const unsigned long &i, const unsigned long &j, const unsigned long &k, const unsigned long &l, const unsigned long &m) const
Range check to catch when an index is out of bounds, if so, it issues a warning message and dies by t...
Definition: matrices.h:2135
unsigned long nindex2() const
Return the range of index 2 of the tensor.
Definition: matrices.h:2492
T & raw_direct_access(const unsigned long &i)
Direct access to internal storage of data in flat-packed C-style column-major format....
Definition: matrices.h:2545
unsigned long nindex5() const
Return the range of index 5 of the tensor.
Definition: matrices.h:2510
unsigned P
3rd Tensor dimension
Definition: matrices.h:2125
unsigned offset(const unsigned long &i, const unsigned long &j, const unsigned long &k) const
Caculate the offset in flat-packed Cy-style, column-major format, required for a given i,...
Definition: matrices.h:2564
void initialise(const T &val)
Initialise all values in the tensor to val.
Definition: matrices.h:2477
unsigned long nindex1() const
Return the range of index 1 of the tensor.
Definition: matrices.h:2486
T * Tensordata
Private internal representation as pointer to data.
Definition: matrices.h:2116
virtual ~RankFiveTensor()
Destructor: delete the pointers.
Definition: matrices.h:2331
unsigned N
1st Tensor dimension
Definition: matrices.h:2119
RankFiveTensor(const unsigned long &n_index1, const unsigned long &n_index2, const unsigned long &n_index3, const unsigned long &n_index4, const unsigned long &n_index5)
Four parameter constructor, general non-square tensor.
Definition: matrices.h:2289
void resize(const unsigned long &n_index1, const unsigned long &n_index2, const unsigned long &n_index3, const unsigned long &n_index4, const unsigned long &n_index5)
Resize to a general tensor.
Definition: matrices.h:2344
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1701
T & raw_direct_access(const unsigned long &i)
Direct access to internal storage of data in flat-packed C-style column-major format....
Definition: matrices.h:2078
unsigned M
2nd Tensor dimension
Definition: matrices.h:1710
unsigned N
1st Tensor dimension
Definition: matrices.h:1707
RankFourTensor()
Empty constructor.
Definition: matrices.h:1769
RankFourTensor(const unsigned long &n_index1, const unsigned long &n_index2, const unsigned long &n_index3, const unsigned long &n_index4)
Four parameter constructor, general non-square tensor.
Definition: matrices.h:1854
void resize(const unsigned long &n_index1, const unsigned long &n_index2, const unsigned long &n_index3, const unsigned long &n_index4, const T &initial_value)
Resize to a general tensor.
Definition: matrices.h:1961
virtual ~RankFourTensor()
Destructor: delete the pointers.
Definition: matrices.h:1892
unsigned offset(const unsigned long &i, const unsigned long &j) const
Caculate the offset in flat-packed C-style, column-major format, required for a given i,...
Definition: matrices.h:2096
T * Tensordata
Private internal representation as pointer to data.
Definition: matrices.h:1704
void resize(const unsigned long &n_index1, const unsigned long &n_index2, const unsigned long &n_index3, const unsigned long &n_index4)
Resize to a general tensor.
Definition: matrices.h:1905
unsigned long nindex4() const
Return the range of index 4 of the tensor.
Definition: matrices.h:2045
unsigned long nindex2() const
Return the range of index 2 of the tensor.
Definition: matrices.h:2033
const T & raw_direct_access(const unsigned long &i) const
Direct access to internal storage of data in flat-packed C-style column-major format....
Definition: matrices.h:2087
unsigned long nindex1() const
Return the range of index 1 of the tensor.
Definition: matrices.h:2027
T & operator()(const unsigned long &i, const unsigned long &j, const unsigned long &k, const unsigned long &l)
Overload the round brackets to give access as a(i,j,k,l)
Definition: matrices.h:2051
T operator()(const unsigned long &i, const unsigned long &j, const unsigned long &k, const unsigned long &l) const
Overload a const version for read-only access as a(i,j,k,l)
Definition: matrices.h:2063
unsigned long nindex3() const
Return the range of index 3 of the tensor.
Definition: matrices.h:2039
RankFourTensor(const unsigned long &n)
One parameter constructor produces a nxnxnxn tensor.
Definition: matrices.h:1838
RankFourTensor(const unsigned long &n_index1, const unsigned long &n_index2, const unsigned long &n_index3, const unsigned long &n_index4, const T &initial_val)
Four parameter constructor, general non-square tensor.
Definition: matrices.h:1874
unsigned P
3rd Tensor dimension
Definition: matrices.h:1713
void range_check(const unsigned long &i, const unsigned long &j, const unsigned long &k, const unsigned long &l) const
Range check to catch when an index is out of bounds, if so, it issues a warning message and dies by t...
Definition: matrices.h:1720
RankFourTensor & operator=(const RankFourTensor &source_tensor)
Copy assignement.
Definition: matrices.h:1801
RankFourTensor(const RankFourTensor &source_tensor)
Copy constructor: Deep copy.
Definition: matrices.h:1772
unsigned Q
4th Tensor dimension
Definition: matrices.h:1716
void initialise(const T &val)
Initialise all values in the tensor to val.
Definition: matrices.h:2018
void resize(const unsigned long &n)
Resize to a square nxnxnxn tensor.
Definition: matrices.h:1899
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1370
virtual ~RankThreeTensor()
Destructor: delete the pointers.
Definition: matrices.h:1532
void resize(const unsigned long &n_index1, const unsigned long &n_index2, const unsigned long &n_index3, const T &initial_value)
Resize to a general tensor.
Definition: matrices.h:1593
unsigned N
1st Tensor dimension
Definition: matrices.h:1376
RankThreeTensor(const unsigned long &n_index1, const unsigned long &n_index2, const unsigned long &n_index3, const T &initial_val)
Three parameter constructor, general non-square tensor.
Definition: matrices.h:1516
unsigned long nindex1() const
Return the range of index 1 of the tensor.
Definition: matrices.h:1651
unsigned long nindex2() const
Return the range of index 2 of the tensor.
Definition: matrices.h:1657
void range_check(const unsigned long &i, const unsigned long &j, const unsigned long &k) const
Range check to catch when an index is out of bounds, if so, it issues a warning message and dies by t...
Definition: matrices.h:1386
RankThreeTensor(const RankThreeTensor &source_tensor)
Copy constructor: Deep copy.
Definition: matrices.h:1428
RankThreeTensor()
Empty constructor.
Definition: matrices.h:1425
unsigned M
2nd Tensor dimension
Definition: matrices.h:1379
T operator()(const unsigned long &i, const unsigned long &j, const unsigned long &k) const
Overload a const version for read-only access as a(i,j,k)
Definition: matrices.h:1680
void resize(const unsigned long &n)
Resize to a square nxnxn tensor.
Definition: matrices.h:1539
void resize(const unsigned long &n_index1, const unsigned long &n_index2, const unsigned long &n_index3)
Resize to a general tensor.
Definition: matrices.h:1545
unsigned long nindex3() const
Return the range of index 3 of the tensor.
Definition: matrices.h:1663
void initialise(const T &val)
Initialise all values in the tensor to val.
Definition: matrices.h:1642
RankThreeTensor & operator=(const RankThreeTensor &source_tensor)
Copy assignement.
Definition: matrices.h:1450
T & operator()(const unsigned long &i, const unsigned long &j, const unsigned long &k)
Overload the round brackets to give access as a(i,j,k)
Definition: matrices.h:1669
RankThreeTensor(const unsigned long &n_index1, const unsigned long &n_index2, const unsigned long &n_index3)
Three parameter constructor, general non-square tensor.
Definition: matrices.h:1498
unsigned P
3rd Tensor dimension
Definition: matrices.h:1382
T * Tensordata
Private internal representation as pointer to data.
Definition: matrices.h:1373
RankThreeTensor(const unsigned long &n)
One parameter constructor produces a cubic nxnxn tensor.
Definition: matrices.h:1483
//////////////////////////////////////////////////////////////// ////////////////////////////////////...
Definition: matrices.h:562
unsigned long Nnz
Number of non-zero values (i.e. size of Value array)
Definition: matrices.h:574
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:628
T * value()
Access to C-style value array.
Definition: matrices.h:616
virtual void output_bottom_right_zero_helper(std::ostream &outfile) const
Output the "bottom right" entry regardless of it being zero or not (this allows automatic detection o...
Definition: matrices.h:648
const T * value() const
Access to C-style value array (const version)
Definition: matrices.h:622
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:634
T * Value
Internal representation of the matrix values, a pointer.
Definition: matrices.h:565
static T Zero
Dummy zero.
Definition: matrices.h:577
virtual void sparse_indexed_output_helper(std::ostream &outfile) const
Indexed output function to print a matrix to the stream outfile as i,j,a(i,j) for a(i,...
Definition: matrices.h:661
void operator=(const SparseMatrix &)=delete
Broken assignment operator.
unsigned long nnz() const
Return the number of nonzero entries.
Definition: matrices.h:640
unsigned long N
Number of rows.
Definition: matrices.h:568
SparseMatrix()
Default constructor.
Definition: matrices.h:581
SparseMatrix(const SparseMatrix &source_matrix)
Copy constructor.
Definition: matrices.h:584
virtual ~SparseMatrix()
Destructor, delete the memory associated with the values.
Definition: matrices.h:609
unsigned long M
Number of columns.
Definition: matrices.h:571
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
double gershgorin_eigenvalue_estimate(const DenseMatrix< CRDoubleMatrix * > &matrix_pt)
Calculates the largest Gershgorin disc whilst preserving the sign. Let A be an n by n matrix,...
Definition: matrices.cc:4003
void deep_copy(const CRDoubleMatrix *const in_matrix_pt, CRDoubleMatrix &out_matrix)
Create a deep copy of the matrix pointed to by in_matrix_pt.
Definition: matrices.h:3490
void concatenate_without_communication(const Vector< LinearAlgebraDistribution * > &row_distribution_pt, const Vector< LinearAlgebraDistribution * > &col_distribution_pt, const DenseMatrix< CRDoubleMatrix * > &matrix_pt, CRDoubleMatrix &result_matrix)
Concatenate CRDoubleMatrix matrices.
Definition: matrices.cc:5223
double inf_norm(const DenseMatrix< CRDoubleMatrix * > &matrix_pt)
Compute infinity (maximum) norm of sub blocks as if it was one matrix.
Definition: matrices.cc:3731
void concatenate(const DenseMatrix< CRDoubleMatrix * > &matrix_pt, CRDoubleMatrix &result_matrix)
Concatenate CRDoubleMatrix matrices. The in matrices are concatenated such that the block structure o...
Definition: matrices.cc:4349
void create_uniformly_distributed_matrix(const unsigned &nrow, const unsigned &ncol, const OomphCommunicator *const comm_pt, const Vector< double > &values, const Vector< int > &column_indices, const Vector< int > &row_start, CRDoubleMatrix &matrix_out)
Builds a uniformly distributed matrix. A locally replicated matrix is constructed then redistributed ...
Definition: matrices.cc:3676
void output()
Doc the command line arguments.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
std::string RayStr
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Create a struct to provide a comparison function for std::sort.
Definition: matrices.h:942
bool operator()(const std::pair< int, double > &pair_1, const std::pair< int, double > &pair_2)
Definition: matrices.h:944