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