oomph::DenseComplexMatrix Class Reference

#include <complex_matrices.h>

Inheritance diagram for oomph::DenseComplexMatrix:

Public Member Functions
	DenseComplexMatrix ()
	Empty constructor, simply assign the lengths N and M to 0. More...
	DenseComplexMatrix (const unsigned long &n)
	Constructor to build a square n by n matrix. More...
	DenseComplexMatrix (const unsigned long &n, const unsigned long &m)
	Constructor to build a matrix with n rows and m columns. More...
	DenseComplexMatrix (const unsigned long &n, const unsigned long &m, const std::complex< double > &initial_val)
	DenseComplexMatrix (const DenseComplexMatrix &matrix)=delete
	Broken copy constructor. More...
void	operator= (const DenseComplexMatrix &)=delete
	Broken assignment operator. More...
unsigned long	nrow () const
	Return the number of rows of the matrix. More...
unsigned long	ncol () const
	Return the number of columns of the matrix. More...
std::complex< double >	operator() (const unsigned long &i, const unsigned long &j) const
std::complex< double > &	operator() (const unsigned long &i, const unsigned long &j)
virtual	~DenseComplexMatrix ()
	Destructor, delete the storage for Index vector and LU storage (if any) More...
int	ludecompose ()
void	lubksub (Vector< std::complex< double >> &rhs)
	Back substitute an LU decomposed matrix. More...
void	residual (const Vector< std::complex< double >> &x, const Vector< std::complex< double >> &rhs, Vector< std::complex< double >> &residual)
	Find the residual of Ax=b, i.e. r=b-Ax. More...
void	multiply (const Vector< std::complex< double >> &x, Vector< std::complex< double >> &soln)
	Multiply the matrix by the vector x: soln=Ax. More...
void	multiply_transpose (const Vector< std::complex< double >> &x, Vector< std::complex< double >> &soln)
	Multiply the transposed matrix by the vector x: soln=A^T x. More...
Public Member Functions inherited from oomph::ComplexMatrixBase
	ComplexMatrixBase ()
	(Empty) constructor. More...
	ComplexMatrixBase (const ComplexMatrixBase &matrix)=delete
	Broken copy constructor. More...
void	operator= (const ComplexMatrixBase &)=delete
	Broken assignment operator. More...
virtual	~ComplexMatrixBase ()
	virtual (empty) destructor More...
virtual void	solve (Vector< std::complex< double >> &rhs)
virtual void	solve (const Vector< std::complex< double >> &rhs, Vector< std::complex< double >> &soln)
virtual double	max_residual (const Vector< std::complex< double >> &x, const Vector< std::complex< double >> &rhs)
Public Member Functions inherited from oomph::DenseMatrix< std::complex< double > >
	DenseMatrix ()
	Empty constructor, simply assign the lengths N and M to 0. More...
	DenseMatrix (const DenseMatrix &source_matrix)
	Copy constructor: Deep copy! More...
	DenseMatrix (const unsigned long &n)
	Constructor to build a square n by n matrix. More...
	DenseMatrix (const unsigned long &n, const unsigned long &m)
	Constructor to build a matrix with n rows and m columns. More...
	DenseMatrix (const unsigned long &n, const unsigned long &m, const std::complex< double > &initial_val)
DenseMatrix &	operator= (const DenseMatrix &source_matrix)
	Copy assignment. More...
std::complex< double > &	entry (const unsigned long &i, const unsigned long &j)
std::complex< double >	get_entry (const unsigned long &i, const unsigned long &j) const
virtual	~DenseMatrix ()
	Destructor, clean up the matrix data. More...
unsigned long	nrow () const
	Return the number of rows of the matrix. More...
unsigned long	ncol () const
	Return the number of columns of the matrix. More...
void	resize (const unsigned long &n)
void	resize (const unsigned long &n, const unsigned long &m)
void	resize (const unsigned long &n, const unsigned long &m, const std::complex< double > &initial_value)
void	initialise (const std::complex< double > &val)
	Initialize all values in the matrix to val. More...
void	output (std::ostream &outfile) const
	Output function to print a matrix row-by-row to the stream outfile. More...
void	output (std::string filename) const
	Output function to print a matrix row-by-row to a file. Specify filename. More...
void	indexed_output (std::ostream &outfile) const
	Indexed output as i,j,a(i,j) More...
void	indexed_output (std::string filename) const
void	output_bottom_right_zero_helper (std::ostream &outfile) const
void	sparse_indexed_output_helper (std::ostream &outfile) const
	Sparse indexed output as i,j,a(i,j) for a(i,j)!=0 only. More...
Public Member Functions inherited from oomph::Matrix< T, MATRIX_TYPE >
	Matrix ()
	(Empty) constructor More...
	Matrix (const Matrix &matrix)=delete
	Broken copy constructor. More...
void	operator= (const Matrix &)=delete
	Broken assignment operator. More...
virtual	~Matrix ()
	Virtual (empty) destructor. More...
T	operator() (const unsigned long &i, const unsigned long &j) const
T &	operator() (const unsigned long &i, const unsigned long &j)
void	sparse_indexed_output (std::ostream &outfile, const unsigned &precision=0, const bool &output_bottom_right_zero=false) const
void	sparse_indexed_output (std::string filename, const unsigned &precision=0, const bool &output_bottom_right_zero=false) const

Private Member Functions
void	delete_lu_factors ()

Private Attributes
Vector< long > *	Index
	Pointer to storage for the index of permutations in the LU solve. More...
std::complex< double > *	LU_factors
	Pointer to storage for the LU decomposition. More...
bool	Overwrite_matrix_storage

Additional Inherited Members
Protected Member Functions inherited from oomph::Matrix< T, MATRIX_TYPE >
void	range_check (const unsigned long &i, const unsigned long &j) const
Protected Attributes inherited from oomph::DenseMatrix< std::complex< double > >
std::complex< double > *	Matrixdata
	Internal representation of matrix as a pointer to data. More...
unsigned long	N
	Number of rows. More...
unsigned long	M
	Number of columns. More...

Detailed Description

Class of matrices containing double complex, and stored as a DenseMatrix<complex<double> >, but with solving functionality inherited from the abstract ComplexMatrix class.

Constructor & Destructor Documentation

DenseComplexMatrix() [1/5]

oomph::DenseComplexMatrix::DenseComplexMatrix ( )

inline

Empty constructor, simply assign the lengths N and M to 0.

      : DenseMatrix<std::complex<double>>(),
        Index(0),
        LU_factors(0),
        Overwrite_matrix_storage(false)
    {
    }

DenseComplexMatrix() [2/5]

oomph::DenseComplexMatrix::DenseComplexMatrix ( const unsigned long &  n )

inline

Constructor to build a square n by n matrix.

      : DenseMatrix<std::complex<double>>(n),
        Index(0),
        LU_factors(0),
        Overwrite_matrix_storage(false)
    {
    }

DenseComplexMatrix() [3/5]

oomph::DenseComplexMatrix::DenseComplexMatrix ( const unsigned long &  n, const unsigned long &  m  )

inline

Constructor to build a matrix with n rows and m columns.

      : DenseMatrix<std::complex<double>>(n, m),
        Index(0),
        LU_factors(0),
        Overwrite_matrix_storage(false)
    {
    }

DenseComplexMatrix() [4/5]

oomph::DenseComplexMatrix::DenseComplexMatrix ( const unsigned long &  n, const unsigned long &  m, const std::complex< double > &  initial_val  )

inline

Constructor to build a matrix with n rows and m columns, with initial value initial_val

      : DenseMatrix<std::complex<double>>(n, m, initial_val),
        Index(0),
        LU_factors(0),
        Overwrite_matrix_storage(false)
    {
    }

DenseComplexMatrix() [5/5]

oomph::DenseComplexMatrix::DenseComplexMatrix ( const DenseComplexMatrix &  matrix )

delete

Broken copy constructor.

~DenseComplexMatrix()

oomph::DenseComplexMatrix::~DenseComplexMatrix ( )

virtual

Destructor, delete the storage for Index vector and LU storage (if any)

Destructor clean up the LU factors that have been allocated.

  {
    // Delete the storage for the index vector
    delete Index;
 
    // Delete the pointer to the LU factors
    delete[] LU_factors;
 
    // Null out the vector
    LU_factors = 0;
  }

References Index, and LU_factors.

Member Function Documentation

delete_lu_factors()

void oomph::DenseComplexMatrix::delete_lu_factors ( )

private

Function that deletes the storage for the LU_factors, if it is not the same as the matrix storage

Delete the storage that has been allocated for the LU factors, if the matrix data is not itself being overwritten.

  {
    // Clean up the LU factor storage, if it has been allocated
    // and it's not the same as the matrix storage
    if ((!Overwrite_matrix_storage) && (LU_factors != 0))
    {
      // Delete the pointer to the LU factors
      delete[] LU_factors;
      // Null out the vector
      LU_factors = 0;
    }
  }

References LU_factors, and Overwrite_matrix_storage.

Referenced by ludecompose().

lubksub()

void oomph::DenseComplexMatrix::lubksub ( Vector< std::complex< double >> &  rhs )

virtual

Back substitute an LU decomposed matrix.

Overload the LU back substitution. Note that the rhs will be overwritten with the solution vector

Reimplemented from oomph::ComplexMatrixBase.

  {
#ifdef PARANOID
    // Check Matrix is square
    if (N != M)
    {
      throw OomphLibError(
        "This matrix is not square,  the matrix MUST be square!",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
    // Check that size of rhs=nrow() and index=nrow() are correct!
    if (rhs.size() != N)
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The rhs vector is not the right size. It is "
                           << rhs.size() << ", it should be " << N << std::endl;
 
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    if (Index == 0)
    {
      throw OomphLibError(
        "Index vector has not been allocated. Have you called ludecompse()\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
    if (Index->size() != N)
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The Index vector is not the right size. It is "
                           << Index->size() << ", it should be " << N
                           << std::endl;
 
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    unsigned long ii = 0;
    for (unsigned i = 0; i < N; i++)
    {
      unsigned long ip = (*Index)[i];
      std::complex<double> sum = rhs[ip];
      rhs[ip] = rhs[i];
 
      if (ii != 0)
      {
        for (unsigned long j = ii - 1; j < i; j++)
        {
          sum -= LU_factors[M * i + j] * rhs[j];
        }
      }
      else if (sum != 0.0)
      {
        ii = i + 1;
      }
      rhs[i] = sum;
    }
 
    // Now do the back substitution
    for (long i = long(N) - 1; i >= 0; i--)
    {
      std::complex<double> sum = rhs[i];
      for (long j = i + 1; j < long(N); j++)
        sum -= LU_factors[M * i + j] * rhs[j];
      rhs[i] = sum / LU_factors[M * i + i];
    }
  }

References i, Index, j, LU_factors, oomph::DenseMatrix< std::complex< double > >::N, OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

ludecompose()

int oomph::DenseComplexMatrix::ludecompose ( )

virtual

Overload the LU decomposition. For this storage scheme the matrix storage will be over-writeen by its LU decomposition

LU decompose a matrix, over-writing it and recording the pivots numbers in the Index vector. Returns the sign of the determinant.

Reimplemented from oomph::ComplexMatrixBase.

  {
#ifdef PARANOID
    // Check Matrix is square
    if (N != M)
    {
      throw OomphLibError(
        "This matrix is not square, the matrix MUST be square!",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Small constant
    const double small_number = 1.0e-20;
 
    // If the Index vector has not already been allocated,
    // allocated storage for the index of permutations
    if (Index == 0)
    {
      Index = new Vector<long>;
    }
 
    // Resize the index vector to the correct length
    Index->resize(N);
 
    // Vector scaling stores the implicit scaling of each row
    Vector<double> scaling(N);
 
    // Integer to store the sign that must multiply the determinant as
    // a consequence of the row/column interchanges
    int signature = 1;
 
    // Loop over rows to get implicit scaling information
    for (unsigned long i = 0; i < N; i++)
    {
      double big = 0.0;
      for (unsigned long j = 0; j < M; j++)
      {
        double tmp = std::abs((*this)(i, j));
        if (tmp > big) big = tmp;
      }
      if (big == 0.0)
      {
        throw OomphLibError(
          "Singular Matrix", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
      // Save the scaling
      scaling[i] = 1.0 / big;
    }
 
    // Firsly, we shall delete any previous LU storage.
    // If the user calls this function twice without changing the matrix
    // then it is their own inefficiency, not ours (this time).
    delete_lu_factors();
 
    // Now if we're saving the LU factors separately (the default).
    if (!Overwrite_matrix_storage)
    {
      // Allocate storage for the LU factors
      // Assign space for the n rows
      LU_factors = new std::complex<double>[N * M];
 
      // Now we know that memory has been allocated, copy over
      // the matrix values
      for (unsigned long i = 0; i < (N * M); i++)
      {
        LU_factors[i] = Matrixdata[i];
      }
    }
    // Otherwise the pointer to the LU_factors is the same as the
    // matrix data
    else
    {
      LU_factors = Matrixdata;
    }
 
    // Loop over columns
    for (unsigned long j = 0; j < M; j++)
    {
      // Initialise imax
      unsigned long imax = 0;
 
      for (unsigned long i = 0; i < j; i++)
      {
        std::complex<double> sum = LU_factors[M * i + j];
        for (unsigned long k = 0; k < i; k++)
        {
          sum -= LU_factors[M * i + k] * LU_factors[M * k + j];
        }
        LU_factors[M * i + j] = sum;
      }
 
      // Initialise search for largest pivot element
      double largest_entry = 0.0;
      for (unsigned long i = j; i < N; i++)
      {
        std::complex<double> sum = LU_factors[M * i + j];
        for (unsigned long k = 0; k < j; k++)
        {
          sum -= LU_factors[M * i + k] * LU_factors[M * k + j];
        }
        LU_factors[M * i + j] = sum;
        // Set temporary
        double tmp = scaling[i] * std::abs(sum);
        if (tmp >= largest_entry)
        {
          largest_entry = tmp;
          imax = i;
        }
      }
 
      // Test to see if we need to interchange rows
      if (j != imax)
      {
        for (unsigned long k = 0; k < N; k++)
        {
          std::complex<double> tmp = LU_factors[M * imax + k];
          LU_factors[M * imax + k] = LU_factors[M * j + k];
          LU_factors[M * j + k] = tmp;
        }
        // Change the parity of signature
        signature = -signature;
        // Interchange scale factor
        scaling[imax] = scaling[j];
      }
 
      // Set the index
      (*Index)[j] = imax;
      if (LU_factors[M * j + j] == 0.0)
      {
        LU_factors[M * j + j] = small_number;
      }
      // Divide by pivot element
      if (j != N - 1)
      {
        std::complex<double> tmp = 1.0 / LU_factors[M * j + j];
        for (unsigned long i = j + 1; i < N; i++)
        {
          LU_factors[M * i + j] *= tmp;
        }
      }
 
    } // End of loop over columns
 
 
    // Now multiply all the diagonal terms together to get the determinant
    // Note that we need to use the mantissa, exponent formulation to
    // avoid underflow errors. This is way more tedious in complex arithmetic.
    std::complex<double> determinant_mantissa(1.0, 0.0);
    std::complex<int> determinant_exponent(0, 0);
    for (unsigned i = 0; i < N; i++)
    {
      // Get the next entry in mantissa exponent form
      std::complex<double> entry;
      int iexp_real = 0, iexp_imag = 0;
      entry =
        std::complex<double>(frexp(LU_factors[M * i + i].real(), &iexp_real),
                             frexp(LU_factors[M * i + i].imag(), &iexp_imag));
 
      // Now calculate the four terms that contribute to the real
      // and imaginary parts
      // i.e. A + Bi  * C + Di = AC - BD + i(AD + BC)
      double temp_mantissa[4];
      int temp_exp[4];
 
      // Get the first contribution to the real part, AC
      temp_mantissa[0] = determinant_mantissa.real() * entry.real();
      temp_exp[0] = determinant_exponent.real() + iexp_real;
      // Get the second contribution to the real part, DB
      temp_mantissa[1] = determinant_mantissa.imag() * entry.imag();
      temp_exp[1] = determinant_exponent.imag() + iexp_imag;
 
      // Get the first contribution to the imaginary part, AD
      temp_mantissa[2] = determinant_mantissa.real() * entry.imag();
      temp_exp[2] = determinant_exponent.real() + iexp_imag;
      // Get the second contribution to the imaginary part, BC
      temp_mantissa[3] = determinant_mantissa.imag() * entry.real();
      temp_exp[3] = determinant_exponent.imag() + iexp_real;
 
      // Align the exponents in the two terms for each sum (real or imaginary)
      // We always align up to the larger exponent
      for (unsigned i = 0; i < 3; i += 2)
      {
        // If the first exponent is smaller, move it up
        if (temp_exp[i] < temp_exp[i + 1])
        {
          // The smaller term
          int lower = temp_exp[i];
          // Loop over the difference in the exponents,
          // scaling down the mantissa each time
          for (int index = lower; index < temp_exp[i + 1]; ++index)
          {
            temp_mantissa[i] /= 2.0;
            ++temp_exp[i];
          }
        }
        // Otherwise the second exponent is smaller
        else
        {
          // The smaller term is the other
          int lower = temp_exp[i + 1];
          // Loop over the difference in the exponents,
          // scaling down the mantissa each time
          for (int index = lower; index < temp_exp[i]; ++index)
          {
            temp_mantissa[i + 1] /= 2.0;
            ++temp_exp[i + 1];
          }
        }
      }
 
      // Now combine the terms AC - BD
      // and Combine the terms AD + BC
      determinant_mantissa =
        std::complex<double>(temp_mantissa[0] - temp_mantissa[1],
                             temp_mantissa[2] + temp_mantissa[3]);
      // The exponents will be the same, so pick one
      determinant_exponent = std::complex<int>(temp_exp[0], temp_exp[2]);
 
      // Finally normalise the result
      determinant_mantissa =
        std::complex<double>(frexp(determinant_mantissa.real(), &iexp_real),
                             frexp(determinant_mantissa.imag(), &iexp_imag));
 
      int temp_real = determinant_exponent.real() + iexp_real;
      int temp_imag = determinant_exponent.imag() + iexp_imag;
 
      determinant_exponent = std::complex<int>(temp_real, temp_imag);
    }
 
    // Integer to store the sign of the determinant
    int sign = 0;
    // Find the sign of the determinant (left or right half-plane)
    if (std::abs(determinant_mantissa) > 0.0)
    {
      sign = 1;
    }
    if (std::abs(determinant_mantissa) < 0.0)
    {
      sign = -1;
    }
 
    // Multiply the sign by the signature
    sign *= signature;
 
    // Return the sign of the determinant
    return sign;
  }

References abs(), delete_lu_factors(), oomph::DenseMatrix< std::complex< double > >::entry(), i, imag(), Index, j, k, helpers::lower(), LU_factors, oomph::DenseMatrix< std::complex< double > >::M, oomph::DenseMatrix< std::complex< double > >::Matrixdata, oomph::DenseMatrix< std::complex< double > >::N, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Overwrite_matrix_storage, SYCL::sign(), and tmp.

multiply()

void oomph::DenseComplexMatrix::multiply ( const Vector< std::complex< double >> &  x, Vector< std::complex< double >> &  soln  )

virtual

Multiply the matrix by the vector x: soln=Ax.

Implements oomph::ComplexMatrixBase.

  {
#ifdef PARANOID
    // Check to see if x.size() = ncol().
    if (x.size() != M)
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The x vector is not the right size. It is "
                           << x.size() << ", it should be " << M << std::endl;
 
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    if (soln.size() != N)
    {
      // Resize and initialize the solution vector
      soln.resize(N);
    }
 
    // Multiply the matrix A, by the vector x
    for (unsigned long i = 0; i < N; i++)
    {
      soln[i] = 0.0;
      for (unsigned long j = 0; j < M; j++)
      {
        soln[i] += Matrixdata[M * i + j] * x[j];
      }
    }
  }

References i, j, oomph::DenseMatrix< std::complex< double > >::M, oomph::DenseMatrix< std::complex< double > >::Matrixdata, oomph::DenseMatrix< std::complex< double > >::N, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and plotDoE::x.

multiply_transpose()

void oomph::DenseComplexMatrix::multiply_transpose ( const Vector< std::complex< double >> &  x, Vector< std::complex< double >> &  soln  )

virtual

Multiply the transposed matrix by the vector x: soln=A^T x.

Implements oomph::ComplexMatrixBase.

  {
#ifdef PARANOID
    // Check to see x.size() = nrow()
    if (x.size() != N)
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The x vector is not the right size. It is "
                           << x.size() << ", it should be " << N << std::endl;
 
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    if (soln.size() != M)
    {
      // Resize and initialize the solution vector
      soln.resize(M);
    }
 
    // Initialise the solution
    for (unsigned long i = 0; i < M; i++)
    {
      soln[i] = 0.0;
    }
 
 
    // Matrix vector product
    for (unsigned long i = 0; i < N; i++)
    {
      for (unsigned long j = 0; j < M; j++)
      {
        soln[j] += Matrixdata[N * i + j] * x[i];
      }
    }
  }

References i, j, oomph::DenseMatrix< std::complex< double > >::M, oomph::DenseMatrix< std::complex< double > >::Matrixdata, oomph::DenseMatrix< std::complex< double > >::N, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and plotDoE::x.

ncol()

unsigned long oomph::DenseComplexMatrix::ncol ( ) const

inlinevirtual

Return the number of columns of the matrix.

Implements oomph::ComplexMatrixBase.

    {
      return DenseMatrix<std::complex<double>>::ncol();
    }

nrow()

unsigned long oomph::DenseComplexMatrix::nrow ( ) const

inlinevirtual

Return the number of rows of the matrix.

Implements oomph::ComplexMatrixBase.

    {
      return DenseMatrix<std::complex<double>>::nrow();
    }

operator()() [1/2]

std::complex<double>& oomph::DenseComplexMatrix::operator() ( const unsigned long &  i, const unsigned long &  j  )

inline

Overload the non-const version of the round-bracket access operator for read-write access

    {
      return DenseMatrix<std::complex<double>>::entry(i, j);
    }

References oomph::DenseMatrix< std::complex< double > >::entry(), i, and j.

operator()() [2/2]

std::complex<double> oomph::DenseComplexMatrix::operator() ( const unsigned long &  i, const unsigned long &  j  ) const

inlinevirtual

Overload the const version of the round-bracket access operator for read-only access.

Implements oomph::ComplexMatrixBase.

    {
      return DenseMatrix<std::complex<double>>::get_entry(i, j);
    }

References oomph::DenseMatrix< std::complex< double > >::get_entry(), i, and j.

operator=()

void oomph::DenseComplexMatrix::operator= ( const DenseComplexMatrix &  )

delete

Broken assignment operator.

residual()

void oomph::DenseComplexMatrix::residual ( const Vector< std::complex< double >> &  x, const Vector< std::complex< double >> &  rhs, Vector< std::complex< double >> &  residual  )

virtual

Find the residual of Ax=b, i.e. r=b-Ax.

Find the residual of Ax=b, ie r=b-Ax for the "solution" x.

Implements oomph::ComplexMatrixBase.

  {
#ifdef PARANOID
    // Check Matrix is square
    if (N != M)
    {
      throw OomphLibError(
        "This matrix is not square, the matrix MUST be square!",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
    // Check that size of rhs = nrow()
    if (rhs.size() != N)
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The rhs vector is not the right size. It is "
                           << rhs.size() << ", it should be " << N << std::endl;
 
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // Check that the size of x is correct
    if (x.size() != N)
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The x vector is not the right size. It is "
                           << x.size() << ", it should be " << N << std::endl;
 
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
    // If size of residual is wrong, correct it!
    if (residual.size() != N)
    {
      residual.resize(N);
    }
 
    // Multiply the matrix by the vector x in residual vector
    for (unsigned long i = 0; i < N; i++)
    {
      residual[i] = rhs[i];
      for (unsigned long j = 0; j < M; j++)
      {
        residual[i] -= Matrixdata[M * i + j] * x[j];
      }
    }
  }

References i, j, oomph::DenseMatrix< std::complex< double > >::M, oomph::DenseMatrix< std::complex< double > >::Matrixdata, oomph::DenseMatrix< std::complex< double > >::N, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and plotDoE::x.

Member Data Documentation

Index

Vector<long>* oomph::DenseComplexMatrix::Index

private

Pointer to storage for the index of permutations in the LU solve.

Referenced by lubksub(), ludecompose(), and ~DenseComplexMatrix().

LU_factors

std::complex<double>* oomph::DenseComplexMatrix::LU_factors

private

Pointer to storage for the LU decomposition.

Referenced by delete_lu_factors(), lubksub(), ludecompose(), and ~DenseComplexMatrix().

Overwrite_matrix_storage

bool oomph::DenseComplexMatrix::Overwrite_matrix_storage

private

Boolean flag used to decided whether the LU decomposition is stored separately, or not

Referenced by delete_lu_factors(), and ludecompose().

---

The documentation for this class was generated from the following files:

• complex_matrices.h
• complex_matrices.cc