DenseMatrix.hpp

Go to the documentation of this file.

/*
 * STRUMPACK -- STRUctured Matrices PACKage, Copyright (c) 2014, The
 * Regents of the University of California, through Lawrence Berkeley
 * National Laboratory (subject to receipt of any required approvals
 * from the U.S. Dept. of Energy).  All rights reserved.
 *
 * If you have questions about your rights to use or distribute this
 * software, please contact Berkeley Lab's Technology Transfer
 * Department at TTD@lbl.gov.
 *
 * NOTICE. This software is owned by the U.S. Department of Energy. As
 * such, the U.S. Government has been granted for itself and others
 * acting on its behalf a paid-up, nonexclusive, irrevocable,
 * worldwide license in the Software to reproduce, prepare derivative
 * works, and perform publicly and display publicly.  Beginning five
 * (5) years after the date permission to assert copyright is obtained
 * from the U.S. Department of Energy, and subject to any subsequent
 * five (5) year renewals, the U.S. Government is granted for itself
 * and others acting on its behalf a paid-up, nonexclusive,
 * irrevocable, worldwide license in the Software to reproduce,
 * prepare derivative works, distribute copies to the public, perform
 * publicly and display publicly, and to permit others to do so.
 *
 * Developers: Pieter Ghysels, Francois-Henry Rouet, Xiaoye S. Li.
 *             (Lawrence Berkeley National Lab, Computational Research
 *             Division).
 *
 */
#ifndef DENSE_MATRIX_HPP
#define DENSE_MATRIX_HPP
 
#include <string>
#include <vector>
#include <functional>
 
#include "misc/RandomWrapper.hpp"
#include "BLASLAPACKWrapper.hpp"
 
 
namespace strumpack {
 
  enum class Trans : char {
    N='N',  
    C='C',  
    T='T'   
  };
 
  inline Trans c2T(char op) {
    switch (op) {
    case 'n': case 'N': return Trans::N;
    case 't': case 'T': return Trans::T;
    case 'c': case 'C': return Trans::C;
    default:
      std::cerr << "ERROR: char " << op << " not recognized,"
                << " should be one of n/N, t/T or c/C" << std::endl;
      return Trans::N;
    }
  }
 
 
  enum class Side : char {
    L='L',  
    R='R'   
  };
 
  enum class UpLo : char {
    U='U',  
    L='L',   
    F='F'   
  };
 
  enum class Diag : char {
    U='U',  
    N='N'   
  };
 
  enum class Jobz : char {
    N='N', 
    V='V'  
  };
 
 
  template<typename scalar_t> class DenseMatrix {
    using real_t = typename RealType<scalar_t>::value_type;
 
  protected:
    scalar_t* data_ = nullptr;
    std::size_t rows_ = 0;
    std::size_t cols_ = 0;
    std::size_t ld_ = 1;
 
  public:
 
    DenseMatrix();
 
    DenseMatrix(std::size_t m, std::size_t n);
 
    DenseMatrix(std::size_t m, std::size_t n,
                const std::function<scalar_t(std::size_t,std::size_t)>& A);
 
    DenseMatrix(std::size_t m, std::size_t n,
                const scalar_t* D, std::size_t ld);
 
    DenseMatrix(std::size_t m, std::size_t n, const DenseMatrix<scalar_t>& D,
                std::size_t i, std::size_t j);
 
    DenseMatrix(const DenseMatrix<scalar_t>& D);
 
    DenseMatrix(DenseMatrix<scalar_t>&& D);
 
    virtual ~DenseMatrix();
 
    virtual DenseMatrix<scalar_t>& operator=(const DenseMatrix<scalar_t>& D);
 
    virtual DenseMatrix<scalar_t>& operator=(DenseMatrix<scalar_t>&& D);
 
    inline std::size_t rows() const { return rows_; }
 
    inline std::size_t cols() const { return cols_; }
 
    inline std::size_t ld() const { return ld_; }
 
    inline const scalar_t* data() const { return data_; }
 
    inline scalar_t* data() { return data_; }
 
    inline scalar_t* end() { return data_ + ld_ * cols_; }
 
    inline const scalar_t* end() const {
      if (rows_ == 0 || cols_ == 0) return data_;
      return data_ + ld_ * cols_;
    }
 
    inline const scalar_t& operator()(std::size_t i, std::size_t j) const
    { assert(i<=rows() && j<=cols()); return data_[i+ld_*j]; }
 
    inline const scalar_t* ptr(std::size_t i, std::size_t j) const
    { assert(i<=rows() && j<=cols()); return data_+i+ld_*j; }
 
    inline scalar_t& operator()(std::size_t i, std::size_t j)
    { assert(i<=rows() && j<=cols()); return data_[i+ld_*j]; }
 
    inline scalar_t* ptr(std::size_t i, std::size_t j)
    { assert(i<=rows() && j<=cols()); return data_+i+ld_*j; }
 
    void print() const { print("A"); }
 
    void print(std::string name, bool all=false, int width=8) const;
 
    void print_to_file(std::string name,
                       std::string filename, int width=8) const;
 
    void random();
 
    void random(random::RandomGeneratorBase<typename RealType<scalar_t>::
                value_type>& rgen);
 
    void fill(scalar_t v);
 
    void fill(const std::function<scalar_t(std::size_t,std::size_t)>& A);
 
    void zero();
 
    void eye();
 
    virtual void clear();
 
    void resize(std::size_t m, std::size_t n);
 
    void hconcat(const DenseMatrix<scalar_t>& b);
 
    void copy(const DenseMatrix<scalar_t>& B,
              std::size_t i=0, std::size_t j=0);
 
    void copy(const scalar_t* B, std::size_t ldb);
 
    DenseMatrix<scalar_t> conj_transpose() const;
 
    void conj_transpose(DenseMatrix<scalar_t>& X) const;
 
    void laswp(const std::vector<int>& P, bool fwd);
 
    void laswp(const int* P, bool fwd);
 
    void lapmr(const std::vector<int>& P, bool fwd);
 
    void lapmt(const std::vector<int>& P, bool fwd);
 
    void extract_rows(const std::vector<std::size_t>& I,
                      DenseMatrix<scalar_t>& B) const;
 
    DenseMatrix<scalar_t>
    extract_rows(const std::vector<std::size_t>& I) const;
 
    void extract_cols(const std::vector<std::size_t>& I,
                      DenseMatrix<scalar_t>& B) const;
 
    DenseMatrix<scalar_t>
    extract_cols(const std::vector<std::size_t>& I) const;
 
    DenseMatrix<scalar_t>
    extract(const std::vector<std::size_t>& I,
            const std::vector<std::size_t>& J) const;
 
    DenseMatrix<scalar_t>&
    scatter_rows_add(const std::vector<std::size_t>& I,
                     const DenseMatrix<scalar_t>& B, int depth);
 
    DenseMatrix<scalar_t>& add(const DenseMatrix<scalar_t>& B, int depth=0);
 
    DenseMatrix<scalar_t>& sub(const DenseMatrix<scalar_t>& B, int depth=0);
 
    DenseMatrix<scalar_t>& scale(scalar_t alpha, int depth=0);
 
    DenseMatrix<scalar_t>&
    scaled_add(scalar_t alpha, const DenseMatrix<scalar_t>& B, int depth=0);
 
    DenseMatrix<scalar_t>&
    scale_and_add(scalar_t alpha, const DenseMatrix<scalar_t>& B, int depth=0);
 
    DenseMatrix<scalar_t>&
    scale_rows(const std::vector<scalar_t>& D, int depth=0);
 
    DenseMatrix<scalar_t>&
    scale_rows_real(const std::vector<real_t>& D, int depth=0);
 
 
    DenseMatrix<scalar_t>& scale_rows(const scalar_t* D, int depth=0);
 
    DenseMatrix<scalar_t>& scale_rows_real(const real_t* D, int depth=0);
 
    DenseMatrix<scalar_t>& div_rows
    (const std::vector<scalar_t>& D, int depth=0);
 
    real_t norm() const;
 
    real_t normF() const;
 
    real_t norm1() const;
 
    real_t normI() const;
 
    virtual std::size_t memory() const {
      return sizeof(scalar_t) * rows() * cols();
    }
 
    virtual std::size_t nonzeros() const {
      return rows()*cols();
    }
 
    int LU(std::vector<int>& piv, int depth=0);
 
    std::vector<int> LU(int depth=0);
 
    int Cholesky(int depth=0);
 
    std::vector<int> LDLt(int depth=0);
 
    //std::vector<int> LDLt_rook(int depth=0);
 
 
    DenseMatrix<scalar_t> solve(const DenseMatrix<scalar_t>& b,
                                const std::vector<int>& piv, int depth=0) const;
 
    void solve_LU_in_place(DenseMatrix<scalar_t>& b,
                           const std::vector<int>& piv, int depth=0) const;
 
    void solve_LU_in_place(DenseMatrix<scalar_t>& b,
                           const int* piv, int depth=0) const;
 
    void solve_LDLt_in_place(DenseMatrix<scalar_t>& b,
                             const std::vector<int>& piv, int depth=0) const;
 
    // void solve_LDLt_rook_in_place
    // (DenseMatrix<scalar_t>& b, const std::vector<int>& piv, int depth=0) const;
 
    void LQ(DenseMatrix<scalar_t>& L,
            DenseMatrix<scalar_t>& Q, int depth) const;
 
    void orthogonalize(scalar_t& r_max, scalar_t& r_min, int depth);
 
    void ID_column(DenseMatrix<scalar_t>& X, std::vector<int>& piv,
                   std::vector<std::size_t>& ind, real_t rel_tol,
                   real_t abs_tol, int max_rank, int depth);
 
    void ID_row(DenseMatrix<scalar_t>& X, std::vector<int>& piv,
                std::vector<std::size_t>& ind, real_t rel_tol, real_t abs_tol,
                int max_rank, int depth) const;
 
    void low_rank(DenseMatrix<scalar_t>& U, DenseMatrix<scalar_t>& V,
                  real_t rel_tol, real_t abs_tol, int max_rank, int depth) const;
 
    std::vector<scalar_t> singular_values() const;
 
    void shift(scalar_t sigma);
 
 
    int syev(Jobz job, UpLo ul, std::vector<scalar_t>& lambda);
 
    void write(const std::string& fname) const;
 
    static DenseMatrix<scalar_t> read(const std::string& fname);
 
 
    std::size_t subnormals() const;
    std::size_t zeros() const;
 
  private:
    void ID_column_GEQP3(DenseMatrix<scalar_t>& X, std::vector<int>& piv,
                         std::vector<std::size_t>& ind, real_t rel_tol,
                         real_t abs_tol, int max_rank, int depth);
 
    template<typename T> friend class DistributedMatrix;
 
    template<typename T> friend std::ofstream&
    operator<<(std::ofstream& os, const DenseMatrix<T>& D);
    template<typename T> friend std::ifstream&
    operator>>(std::ifstream& is, DenseMatrix<T>& D);
  };
 
 
  template<typename scalar_t>
  class DenseMatrixWrapper : public DenseMatrix<scalar_t> {
  public:
    DenseMatrixWrapper() : DenseMatrix<scalar_t>() {}
 
    DenseMatrixWrapper(std::size_t m, std::size_t n,
                       scalar_t* D, std::size_t ld) noexcept {
      this->data_ = D; this->rows_ = m; this->cols_ = n;
      this->ld_ = std::max(std::size_t(1), ld);
    }
 
    DenseMatrixWrapper(std::size_t m, std::size_t n, DenseMatrix<scalar_t>& D,
                       std::size_t i, std::size_t j) noexcept
      : DenseMatrixWrapper<scalar_t>(m, n, &D(i, j), D.ld()) {
      assert(i+m <= D.rows());
      assert(j+n <= D.cols());
    }
 
    DenseMatrixWrapper(DenseMatrixWrapper<scalar_t>& D) noexcept {
      this->data_ = D.data();
      this->rows_ = D.rows();
      this->cols_ = D.cols();
      this->ld_ = D.ld();
    }
 
    DenseMatrixWrapper(DenseMatrixWrapper<scalar_t>&& D) noexcept {
      this->data_ = D.data();
      this->rows_ = D.rows();
      this->cols_ = D.cols();
      this->ld_ = D.ld();
    }
 
    DenseMatrixWrapper(const DenseMatrix<scalar_t>&) = delete;
    DenseMatrixWrapper(const DenseMatrixWrapper<scalar_t>& D) = delete;
 
    DenseMatrixWrapper(DenseMatrix<scalar_t>&&) = delete;
 
    virtual ~DenseMatrixWrapper() { this->data_ = nullptr; }
 
    void clear() override {
      this->rows_ = 0; this->cols_ = 0;
      this->ld_ = 1; this->data_ = nullptr;
    }
 
    std::size_t memory() const override { return 0; }
 
    std::size_t nonzeros() const override { return 0; }
 
 
    DenseMatrixWrapper<scalar_t>&
    operator=(const DenseMatrixWrapper<scalar_t>& D) = delete;
 
    DenseMatrixWrapper<scalar_t>&
    operator=(DenseMatrixWrapper<scalar_t>&& D) {
      this->data_ = D.data(); this->rows_ = D.rows();
      this->cols_ = D.cols(); this->ld_ = D.ld(); return *this; }
 
    DenseMatrix<scalar_t>&
    operator=(const DenseMatrix<scalar_t>& a) override {
      assert(a.rows()==this->rows() && a.cols()==this->cols());
      for (std::size_t j=0; j<this->cols(); j++)
        for (std::size_t i=0; i<this->rows(); i++)
          this->operator()(i, j) = a(i, j);
      return *this;
    }
 
    using DenseMatrix<scalar_t>::operator=;
  };
 
 
  template<typename scalar_t>
  std::unique_ptr<const DenseMatrixWrapper<scalar_t>>
  ConstDenseMatrixWrapperPtr
  (std::size_t m, std::size_t n, const scalar_t* D, std::size_t ld) {
    return std::unique_ptr<const DenseMatrixWrapper<scalar_t>>
      (new DenseMatrixWrapper<scalar_t>(m, n, const_cast<scalar_t*>(D), ld));
  }
 
  template<typename scalar_t>
  std::unique_ptr<const DenseMatrixWrapper<scalar_t>>
  ConstDenseMatrixWrapperPtr
  (std::size_t m, std::size_t n, const DenseMatrix<scalar_t>& D,
   std::size_t i, std::size_t j) {
    return std::unique_ptr<const DenseMatrixWrapper<scalar_t>>
      (new DenseMatrixWrapper<scalar_t>
       (m, n, const_cast<DenseMatrix<scalar_t>&>(D), i, j));
  }
 
 
  template<typename scalar_from_t, typename scalar_to_t> void
  copy(std::size_t m, std::size_t n, const DenseMatrix<scalar_from_t>& a,
       std::size_t ia, std::size_t ja, DenseMatrix<scalar_to_t>& b,
       std::size_t ib, std::size_t jb) {
    for (std::size_t j=0; j<n; j++)
      for (std::size_t i=0; i<m; i++)
        b(ib+i, jb+j) = static_cast<scalar_to_t>(a(ia+i, ja+j));
  }
 
  template<typename scalar_from_t, typename scalar_to_t> void
  copy(const DenseMatrix<scalar_from_t>& a, DenseMatrix<scalar_to_t>& b,
       std::size_t ib=0, std::size_t jb=0) {
    copy(a.rows(), a.cols(), a, 0, 0, b, ib, jb);
  }
 
  template<typename scalar_t> void
  copy(const DenseMatrix<scalar_t>& a, scalar_t* b, std::size_t ldb) {
    for (std::size_t j=0; j<a.cols(); j++)
      for (std::size_t i=0; i<a.rows(); i++)
        b[i+j*ldb] = a(i, j);
  }
 
 
  template<typename scalar_t> DenseMatrix<scalar_t>
  vconcat(const DenseMatrix<scalar_t>& a, const DenseMatrix<scalar_t>& b) {
    assert(a.cols() == b.cols());
    DenseMatrix<scalar_t> tmp(a.rows()+b.rows(), a.cols());
    copy(a, tmp, 0, 0);
    copy(b, tmp, a.rows(), 0);
    return tmp;
  }
 
  template<typename scalar_t> DenseMatrix<scalar_t>
  hconcat(const DenseMatrix<scalar_t>& a, const DenseMatrix<scalar_t>& b) {
    assert(a.rows() == b.rows());
    DenseMatrix<scalar_t> tmp(a.rows(), a.cols()+b.cols());
    copy(a, tmp, 0, 0);
    copy(b, tmp, 0, a.cols());
    return tmp;
  }
 
  template<typename scalar_t> DenseMatrix<scalar_t>
  eye(std::size_t m, std::size_t n) {
    DenseMatrix<scalar_t> I(m, n);
    I.eye();
    return I;
  }
 
 
 
 
  template<typename scalar_t> void
  gemm(Trans ta, Trans tb, scalar_t alpha, const DenseMatrix<scalar_t>& a,
       const DenseMatrix<scalar_t>& b, scalar_t beta,
       DenseMatrix<scalar_t>& c, int depth=0);
 
  template<typename scalar_t> void
  gemm(Trans ta, Trans tb, scalar_t alpha, const DenseMatrix<scalar_t>& a,
       const scalar_t* b, int ldb, scalar_t beta,
       DenseMatrix<scalar_t>& c, int depth=0);
 
  template<typename scalar_t> void
  gemm(Trans ta, Trans tb, scalar_t alpha, const DenseMatrix<scalar_t>& a,
       const DenseMatrix<scalar_t>& b, scalar_t beta,
       scalar_t* c, int ldc, int depth=0);
 
  template<typename scalar_t> void
  trmm(Side s, UpLo ul, Trans ta, Diag d, scalar_t alpha,
       const DenseMatrix<scalar_t>& a, DenseMatrix<scalar_t>& b,
       int depth=0);
 
  template<typename scalar_t> void
  trsm(Side s, UpLo ul, Trans ta, Diag d, scalar_t alpha,
       const DenseMatrix<scalar_t>& a, DenseMatrix<scalar_t>& b,
       int depth=0);
 
  template<typename scalar_t> void
  trsv(UpLo ul, Trans ta, Diag d, const DenseMatrix<scalar_t>& a,
       DenseMatrix<scalar_t>& b, int depth=0);
 
  template<typename scalar_t> void
  gemv(Trans ta, scalar_t alpha, const DenseMatrix<scalar_t>& a,
       const DenseMatrix<scalar_t>& x, scalar_t beta,
       DenseMatrix<scalar_t>& y, int depth=0);
 
  template<typename scalar_t> void
  gemv(Trans ta, scalar_t alpha, const DenseMatrix<scalar_t>& a,
       const scalar_t* x, int incx, scalar_t beta,
       DenseMatrix<scalar_t>& y, int depth=0);
 
  template<typename scalar_t> void
  gemv(Trans ta, scalar_t alpha, const DenseMatrix<scalar_t>& a,
       const DenseMatrix<scalar_t>& x, scalar_t beta,
       scalar_t* y, int incy, int depth=0);
 
  template<typename scalar_t> void
  gemv(Trans ta, scalar_t alpha, const DenseMatrix<scalar_t>& a,
       const scalar_t* x, int incx, scalar_t beta,
       scalar_t* y, int incy, int depth=0);
 
 
  template<typename scalar_t> long long int
  LU_flops(const DenseMatrix<scalar_t>& a) {
    return (is_complex<scalar_t>() ? 4:1) *
      blas::getrf_flops(a.rows(), a.cols());
  }
 
  template<typename scalar_t> long long int
  solve_flops(const DenseMatrix<scalar_t>& b) {
    return (is_complex<scalar_t>() ? 4:1) *
      blas::getrs_flops(b.rows(), b.cols());
  }
 
  template<typename scalar_t> long long int
  LQ_flops(const DenseMatrix<scalar_t>& a) {
    auto minrc = std::min(a.rows(), a.cols());
    return (is_complex<scalar_t>() ? 4:1) *
      (blas::gelqf_flops(a.rows(), a.cols()) +
       blas::xxglq_flops(a.cols(), a.cols(), minrc));
  }
 
  template<typename scalar_t> long long int
  ID_row_flops(const DenseMatrix<scalar_t>& a, int rank) {
    return (is_complex<scalar_t>() ? 4:1) *
      (blas::geqp3_flops(a.cols(), a.rows()) +
       blas::trsm_flops(rank, a.cols() - rank, scalar_t(1.), 'L'));
  }
 
  template<typename scalar_t> long long int
  trsm_flops(Side s, scalar_t alpha, const DenseMatrix<scalar_t>& a,
             const DenseMatrix<scalar_t>& b) {
    return (is_complex<scalar_t>() ? 4:1) *
      blas::trsm_flops(b.rows(), b.cols(), alpha, char(s));
  }
 
  template<typename scalar_t> long long int
  gemm_flops(Trans ta, Trans tb, scalar_t alpha,
             const DenseMatrix<scalar_t>& a,
             const DenseMatrix<scalar_t>& b, scalar_t beta) {
    return (is_complex<scalar_t>() ? 4:1) *
      blas::gemm_flops
      ((ta==Trans::N) ? a.rows() : a.cols(),
       (tb==Trans::N) ? b.cols() : b.rows(),
       (ta==Trans::N) ? a.cols() : a.rows(), alpha, beta);
  }
 
  template<typename scalar_t> long long int
  gemm_flops(Trans ta, Trans tb, scalar_t alpha,
             const DenseMatrix<scalar_t>& a, scalar_t beta,
             const DenseMatrix<scalar_t>& c) {
    return (is_complex<scalar_t>() ? 4:1) *
      blas::gemm_flops
      (c.rows(), c.cols(), (ta==Trans::N) ? a.cols() : a.rows(), alpha, beta);
  }
 
  template<typename scalar_t> long long int
  orthogonalize_flops(const DenseMatrix<scalar_t>& a) {
    auto minrc = std::min(a.rows(), a.cols());
    return (is_complex<scalar_t>() ? 4:1) *
      (blas::geqrf_flops(a.rows(), minrc) +
       blas::xxgqr_flops(a.rows(), minrc, minrc));
  }
 
  template<typename scalar_t,typename cast_t>
  DenseMatrix<cast_t> cast_matrix(const DenseMatrix<scalar_t>& mat);
 
} // end namespace strumpack
 
#endif // DENSE_MATRIX_HPP