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