complex_matrices.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// This header file contains classes and inline function definitions for
// matrices of complex numbers and their derived types
 
// Include guards to prevent multiple inclusion of the header
#ifndef OOMPH_COMPLEX_MATRICES_HEADER
#define OOMPH_COMPLEX_MATRICES_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif
 
// Of course we need the complex type
#include <complex>
 
// Also need the standard matrices header
#include "matrices.h"
 
namespace oomph
{
  //===================================================================
  //===================================================================
  class ComplexMatrixBase
  {
  public:
    ComplexMatrixBase() {}
 
    ComplexMatrixBase(const ComplexMatrixBase& matrix) = delete;
 
    void operator=(const ComplexMatrixBase&) = delete;
 
    virtual unsigned long nrow() const = 0;
 
    virtual unsigned long ncol() const = 0;
 
    virtual ~ComplexMatrixBase() {}
 
    virtual std::complex<double> operator()(const unsigned long& i,
                                            const unsigned long& j) const = 0;
 
    virtual int ludecompose()
    {
      throw OomphLibError(
        "ludecompose() has not been written for this matrix class\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
 
      return 1;
    }
 
    virtual void lubksub(Vector<std::complex<double>>& rhs)
    {
      throw OomphLibError(
        "lubksub() has not been written for this matrix class\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
 
    virtual void solve(Vector<std::complex<double>>& rhs);
 
    virtual void solve(const Vector<std::complex<double>>& rhs,
                       Vector<std::complex<double>>& soln);
 
    virtual void residual(const Vector<std::complex<double>>& x,
                          const Vector<std::complex<double>>& b,
                          Vector<std::complex<double>>& residual) = 0;
 
    virtual double max_residual(const Vector<std::complex<double>>& x,
                                const Vector<std::complex<double>>& rhs)
    {
      unsigned long n = rhs.size();
      Vector<std::complex<double>> res(n);
      residual(x, rhs, res);
      double ans = 0.0;
      for (unsigned long i = 0; i < n; i++)
      {
        ans = std::max(ans, std::abs(res[i]));
      }
      return ans;
    }
 
    virtual void multiply(const Vector<std::complex<double>>& x,
                          Vector<std::complex<double>>& soln) = 0;
 
    virtual void multiply_transpose(const Vector<std::complex<double>>& x,
                                    Vector<std::complex<double>>& soln) = 0;
  };
 
 
  //=================================================================
  //=================================================================
  class DenseComplexMatrix : public ComplexMatrixBase,
                             public DenseMatrix<std::complex<double>>
  {
  private:
    Vector<long>* Index;
 
    std::complex<double>* LU_factors;
 
    bool Overwrite_matrix_storage;
 
    void delete_lu_factors();
 
  public:
    DenseComplexMatrix()
      : DenseMatrix<std::complex<double>>(),
        Index(0),
        LU_factors(0),
        Overwrite_matrix_storage(false)
    {
    }
 
    DenseComplexMatrix(const unsigned long& n)
      : DenseMatrix<std::complex<double>>(n),
        Index(0),
        LU_factors(0),
        Overwrite_matrix_storage(false)
    {
    }
 
    DenseComplexMatrix(const unsigned long& n, const unsigned long& m)
      : DenseMatrix<std::complex<double>>(n, m),
        Index(0),
        LU_factors(0),
        Overwrite_matrix_storage(false)
    {
    }
 
    DenseComplexMatrix(const unsigned long& n,
                       const unsigned long& m,
                       const std::complex<double>& initial_val)
      : DenseMatrix<std::complex<double>>(n, m, initial_val),
        Index(0),
        LU_factors(0),
        Overwrite_matrix_storage(false)
    {
    }
 
 
    DenseComplexMatrix(const DenseComplexMatrix& matrix) = delete;
 
    void operator=(const DenseComplexMatrix&) = delete;
 
    inline unsigned long nrow() const
    {
      return DenseMatrix<std::complex<double>>::nrow();
    }
 
    inline unsigned long ncol() const
    {
      return DenseMatrix<std::complex<double>>::ncol();
    }
 
    inline std::complex<double> operator()(const unsigned long& i,
                                           const unsigned long& j) const
    {
      return DenseMatrix<std::complex<double>>::get_entry(i, j);
    }
 
    inline std::complex<double>& operator()(const unsigned long& i,
                                            const unsigned long& j)
    {
      return DenseMatrix<std::complex<double>>::entry(i, j);
    }
 
    virtual ~DenseComplexMatrix();
 
    int ludecompose();
 
    void lubksub(Vector<std::complex<double>>& rhs);
 
    void residual(const Vector<std::complex<double>>& x,
                  const Vector<std::complex<double>>& rhs,
                  Vector<std::complex<double>>& residual);
 
    void multiply(const Vector<std::complex<double>>& x,
                  Vector<std::complex<double>>& soln);
 
    void multiply_transpose(const Vector<std::complex<double>>& x,
                            Vector<std::complex<double>>& soln);
  };
 
 
  //=================================================================
  //=================================================================
  class CRComplexMatrix : public CRMatrix<std::complex<double>>,
                          public ComplexMatrixBase
  {
  private:
    void* F_factors;
 
    int Info;
 
  public:
    CRComplexMatrix() : CRMatrix<std::complex<double>>(), F_factors(0), Info(0)
    {
      // By default SuperLU solves linear systems quietly
      Doc_stats_during_solve = false;
    }
 
    CRComplexMatrix(const Vector<std::complex<double>>& value,
                    const Vector<int>& column_index,
                    const Vector<int>& row_start,
                    const unsigned long& n,
                    const unsigned long& m)
      : CRMatrix<std::complex<double>>(value, column_index, row_start, n, m),
        F_factors(0),
        Info(0)
    {
      // By default SuperLU solves linear systems quietly
      Doc_stats_during_solve = false;
    }
 
 
    CRComplexMatrix(const CRComplexMatrix& matrix) = delete;
 
    void operator=(const CRComplexMatrix&) = delete;
 
    virtual ~CRComplexMatrix()
    {
      clean_up_memory();
    }
 
    void enable_doc_stats()
    {
      Doc_stats_during_solve = true;
    }
 
    // Set flag to indicate that stats are not to be displayed during
    void disable_doc_stats()
    {
      Doc_stats_during_solve = false;
    }
 
    inline unsigned long nrow() const
    {
      return CRMatrix<std::complex<double>>::nrow();
    }
 
    inline unsigned long ncol() const
    {
      return CRMatrix<std::complex<double>>::ncol();
    }
 
    inline std::complex<double> operator()(const unsigned long& i,
                                           const unsigned long& j) const
    {
      return CRMatrix<std::complex<double>>::get_entry(i, j);
    }
 
    int ludecompose();
 
    void lubksub(Vector<std::complex<double>>& rhs);
 
    void clean_up_memory();
 
    void residual(const Vector<std::complex<double>>& x,
                  const Vector<std::complex<double>>& b,
                  Vector<std::complex<double>>& residual);
 
    void multiply(const Vector<std::complex<double>>& x,
                  Vector<std::complex<double>>& soln);
 
 
    void multiply_transpose(const Vector<std::complex<double>>& x,
                            Vector<std::complex<double>>& soln);
 
  protected:
    bool Doc_stats_during_solve;
  };
 
 
  //=================================================================
  //=================================================================
  class CCComplexMatrix : public ComplexMatrixBase,
                          public CCMatrix<std::complex<double>>
  {
  private:
    void* F_factors;
 
    int Info;
 
  protected:
    bool Doc_stats_during_solve;
 
  public:
    CCComplexMatrix() : CCMatrix<std::complex<double>>(), F_factors(0), Info(0)
    {
      // By default SuperLU solves linear systems quietly
      Doc_stats_during_solve = false;
    }
 
    CCComplexMatrix(const Vector<std::complex<double>>& value,
                    const Vector<int>& row_index,
                    const Vector<int>& column_start,
                    const unsigned long& n,
                    const unsigned long& m)
      : CCMatrix<std::complex<double>>(value, row_index, column_start, n, m),
        F_factors(0),
        Info(0)
    {
      // By default SuperLU solves linear systems quietly
      Doc_stats_during_solve = false;
    }
 
    CCComplexMatrix(const CCComplexMatrix& matrix) = delete;
 
    void operator=(const CCComplexMatrix&) = delete;
 
    virtual ~CCComplexMatrix()
    {
      clean_up_memory();
    }
 
    void enable_doc_stats()
    {
      Doc_stats_during_solve = true;
    }
 
    // Set flag to indicate that stats are not to be displayed during
    void disable_doc_stats()
    {
      Doc_stats_during_solve = false;
    }
 
    inline unsigned long nrow() const
    {
      return CCMatrix<std::complex<double>>::nrow();
    }
 
    inline unsigned long ncol() const
    {
      return CCMatrix<std::complex<double>>::ncol();
    }
 
    inline std::complex<double> operator()(const unsigned long& i,
                                           const unsigned long& j) const
    {
      return CCMatrix<std::complex<double>>::get_entry(i, j);
    }
 
    int ludecompose();
 
    void lubksub(Vector<std::complex<double>>& rhs);
 
    void clean_up_memory();
 
    void residual(const Vector<std::complex<double>>& x,
                  const Vector<std::complex<double>>& b,
                  Vector<std::complex<double>>& residual);
 
 
    void multiply(const Vector<std::complex<double>>& x,
                  Vector<std::complex<double>>& soln);
 
 
    void multiply_transpose(const Vector<std::complex<double>>& x,
                            Vector<std::complex<double>>& soln);
  };
} // namespace oomph
 
#endif