al-matrix.H

Go to the documentation of this file.

/*
                          Aleph_w
 
  Data structures & Algorithms
  version 2.0.0b
  https://github.com/lrleon/Aleph-w
 
  This file is part of Aleph-w library
 
  Copyright (c) 2002-2026 Leandro Rabindranath Leon
 
  Permission is hereby granted, free of charge, to any person obtaining a copy
  of this software and associated documentation files (the "Software"), to deal
  in the Software without restriction, including without limitation the rights
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  copies of the Software, and to permit persons to whom the Software is
  furnished to do so, subject to the following conditions:
 
  The above copyright notice and this permission notice shall be included in all
  copies or substantial portions of the Software.
 
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  SOFTWARE.
*/
 
 
# ifndef AH_MATRIX_H
# define AH_MATRIX_H
 
# include <al-vector.H>
# include <ah-errors.H>
 
namespace Aleph
{
  template <typename Trow = int, typename Tcol = int, typename NumType = double>
  class Matrix
  {
  public:
    using RDomain = AlDomain<Trow>;
    using CDomain = AlDomain<Tcol>;
    using RDomainPtr = std::shared_ptr<RDomain>;
    using CDomainPtr = std::shared_ptr<CDomain>;
 
    using RCPair = std::pair<Trow, Tcol>;
 
    using Map = Aleph::HashMap<RCPair, NumType, MapODhash>;
 
    using ColVector = Vector<Trow, NumType>;
 
    using RowVector = Vector<Tcol, NumType>;
 
  private:
    static const NumType default_epsilon;
 
    RDomainPtr rdomain_ptr = nullptr;
    CDomainPtr cdomain_ptr = nullptr;
    Map entries;
    NumType epsilon;
 
    bool is_zero(const NumType & val) const noexcept
    {
      assert(epsilon >= 0);
      return val >= -epsilon and val <= epsilon;
    }
 
  public:
    const NumType &get_epsilon() const noexcept { return epsilon; }
 
    void set_epsilon(const NumType & e)
    {
      ah_domain_error_if(e < 0)
        << "epsilon must be greater or equal to zero";
      epsilon = e;
    }
 
    typedef std::pair<RCPair, NumType> Pair;
 
    const RDomain &get_row_domain() const noexcept { return *rdomain_ptr; }
 
    const CDomain &get_col_domain() const noexcept { return *cdomain_ptr; }
 
    const RDomainPtr &get_row_domain_ptr() const noexcept { return rdomain_ptr; }
 
    const CDomainPtr &get_col_domain_ptr() const noexcept { return cdomain_ptr; }
 
    Matrix(RDomainPtr rdomain, CDomainPtr cdomain,
           const NumType & e = default_epsilon)
      : rdomain_ptr(std::move(rdomain)),
        cdomain_ptr(std::move(cdomain)),
        entries(10, Aleph::pair_dft_hash_fct<Trow, Tcol>,
                Aleph::pair_snd_hash_fct<Trow, Tcol>)
    {
      set_epsilon(e);
    }
 
    Matrix(const RDomain & rdomain, const CDomain & cdomain,
           const NumType & e = default_epsilon)
      : rdomain_ptr(RDomainPtr(&const_cast<RDomain &>(rdomain), [](RDomain *) {})),
        cdomain_ptr(CDomainPtr(&const_cast<CDomain &>(cdomain), [](CDomain *) {})),
        entries(10, Aleph::pair_dft_hash_fct<Trow, Tcol>,
                Aleph::pair_snd_hash_fct<Trow, Tcol>)
    {
      set_epsilon(e);
    }
 
    Matrix(RDomainPtr rdomain, CDomainPtr cdomain,
           std::initializer_list<std::initializer_list<NumType>> l,
           const NumType & e = default_epsilon)
      : Matrix(std::move(rdomain), std::move(cdomain), e)
    {
      ah_range_error_if(l.size() != rdomain_ptr->size()) << "mismatch with number of rows";
 
      auto itr = l.begin();
      const DynList<Tcol> clist = sort(cdomain_ptr->keys());
      sort(rdomain_ptr->keys()).
          for_each([this, &itr, &clist](const Trow & row)
                     {
                       ah_range_error_if(itr->size() != cdomain_ptr->size())
                          << "mismatch with number of columns";
                       auto itc = itr->begin();
                       clist.for_each([this, &itc, &row](const Tcol & col)
                                        {
                                          const NumType & val = *itc++;
                                          if (is_zero(val))
                                            return;
 
                                          entries.insert(RCPair(row, col), val);
                                        });
                       ++itr;
                     });
    }
 
    Matrix(const RDomain & rdomain, const CDomain & cdomain,
           std::initializer_list<std::initializer_list<NumType>> l,
           const NumType & e = default_epsilon)
      : Matrix(RDomainPtr(&const_cast<RDomain &>(rdomain), [](RDomain *) {}),
               CDomainPtr(&const_cast<CDomain &>(cdomain), [](CDomain *) {}), l, e)
    {}
 
    Matrix(RDomainPtr rdomain, CDomainPtr cdomain,
           const DynList<DynList<NumType>> & l,
           const NumType & e = default_epsilon)
      : Matrix(std::move(rdomain), std::move(cdomain), e)
    {
      ah_range_error_if(l.size() != rdomain_ptr->size()) << "mismatch with number of rows";
 
      for (auto it = zip_it(rdomain_ptr->keys(), l); it.has_curr(); it.next())
        {
          auto t = it.get_curr();
          const Trow & row_idx = get<0>(t);
          const DynList<NumType> & row = get<1>(t);
          for (auto it = zip_it(cdomain_ptr->keys(), row); it.has_curr(); it.next())
            {
              auto t = it.get_curr();
              const NumType & val = get<1>(t);
              if (is_zero(val))
                continue;
              const Tcol & col_idx = get<0>(t);
              entries.insert(RCPair(row_idx, col_idx), val);
            }
        }
    }
 
    Matrix(const RDomain & rdomain, const CDomain & cdomain,
           const DynList<DynList<NumType>> & l,
           const NumType & e = default_epsilon)
      : Matrix(RDomainPtr(&const_cast<RDomain &>(rdomain), [](RDomain *) {}),
               CDomainPtr(&const_cast<CDomain &>(cdomain), [](CDomain *) {}), l, e)
    {}
 
    Matrix(const Matrix & other)
      : rdomain_ptr(other.rdomain_ptr),
        cdomain_ptr(other.cdomain_ptr),
        entries(other.entries),
        epsilon(other.epsilon)
    {
      // empty
    }
 
    Matrix(Matrix && other) noexcept
      : rdomain_ptr(other.rdomain_ptr),
        cdomain_ptr(other.cdomain_ptr),
        entries(std::move(other.entries)),
        epsilon(other.epsilon)
    {
      // empty
    }
 
    Matrix &operator =(const Matrix & other)
    {
      if (this != &other)
        {
          rdomain_ptr = other.rdomain_ptr;
          cdomain_ptr = other.cdomain_ptr;
          entries = other.entries;
          epsilon = other.epsilon;
        }
      return *this;
    }
 
    Matrix &operator =(Matrix && other) noexcept
    {
      if (this != &other)
        {
          rdomain_ptr = other.rdomain_ptr;
          cdomain_ptr = other.cdomain_ptr;
          entries = std::move(other.entries);
          epsilon = other.epsilon;
        }
      return *this;
    }
 
    static Matrix create_by_rows(const RDomain & rdomain,
                                 const DynList<Vector<Tcol, NumType>> & rows,
                                 const NumType & e = default_epsilon)
    {
      Matrix ret(rdomain, rows.get_first().get_domain(), e);
 
      for (auto it = zip_it(rdomain.keys(), rows); it.has_curr(); it.next())
        {
          auto t = it.get_curr();
          const Trow & row_idx = get<0>(t);
          const Vector<Tcol, NumType> & row = get<1>(t);
          for (auto it_col = row.get_it(); it_col.has_curr(); it_col.next())
            {
              auto & p = it_col.get_curr();
              const NumType & val = p.second;
              if (ret.is_zero(val))
                continue;
              const Tcol & col_idx = p.first;
              ret.entries.insert(RCPair(row_idx, col_idx), val);
            }
        }
 
      return ret;
    }
 
    static Matrix create_by_columns(const CDomain & cdomain,
                                    const DynList<Vector<Trow, NumType>> & cols,
                                    const NumType & e = default_epsilon)
    {
      Matrix ret(cols.get_first().get_domain(), cdomain, e);
 
      for (auto it = zip_it(cdomain.keys(), cols); it.has_curr(); it.next())
        {
          auto t = it.get_curr();
          const Tcol & col_idx = get<0>(t);
          const Vector<Trow, NumType> & col = get<1>(t);
          for (auto it_col = col.get_it(); it_col.has_curr(); it_col.next())
            {
              auto & p = it_col.get_curr();
              const NumType & val = p.second;
              if (ret.is_zero(val))
                continue;
              const Trow & row_idx = p.first;
              ret.entries.insert(RCPair(row_idx, col_idx), val);
            }
        }
 
      return ret;
    }
 
    NumType get_entry(const Trow & row, const Tcol & col)
    {
      assert(rdomain_ptr->has(row) and cdomain_ptr->has(col));
      auto *ptr = entries.search(RCPair(row, col));
      if (ptr == nullptr)
        return 0;
      if (is_zero(ptr->second))
        {
          entries.remove_ptr(ptr);
          return 0;
        }
      return ptr->second;
    }
 
    NumType get_entry(const Trow & row, const Tcol & col) const noexcept
    {
      assert(rdomain_ptr->has(row) and cdomain_ptr->has(col));
      auto *ptr = entries.search(RCPair(row, col));
      if (ptr == nullptr)
        return 0;
      return ptr->second;
    }
 
    void set_entry(const Trow & row, const Tcol & col, const NumType & val)
    {
      assert(rdomain_ptr->has(row) and cdomain_ptr->has(col));
      auto *ptr = entries.search(RCPair(row, col));
      if (is_zero(val))
        {
          if (ptr != nullptr)
            entries.remove_ptr(ptr);
          return;
        }
      if (ptr == nullptr)
        entries.insert(RCPair(row, col), val);
      else
        ptr->second = val;
    }
 
    Matrix transpose() const
    {
      Matrix ret_val(*cdomain_ptr, *rdomain_ptr);
      entries.for_each([&ret_val](const Pair & p)
                         {
                           const RCPair & rc = p.first;
                           const NumType & val = p.second;
                           ret_val.set_entry(rc.second, rc.first, val);
                         });
      return ret_val;
    }
 
    DynList<Trow> row_domain_list() const
    {
      return sort(rdomain_ptr->keys());
    }
 
    DynList<Tcol> col_domain_list() const
    {
      return sort(cdomain_ptr->keys());
    }
 
    DynList<NumType> get_row_as_list(const Trow & row) const
    {
      ah_domain_error_if(not rdomain_ptr->has(row)) << "row is not inside the row domain";
      return col_domain_list().
          template maps<NumType>([this, &row](const Tcol & col)
                                   {
                                     return get_entry(row, col);
                                   });
    }
 
    DynList<NumType> get_col_as_list(const Tcol & col) const
    {
      ah_domain_error_if(not cdomain_ptr->has(col)) << "col is not inside the column domain";
      return row_domain_list().template maps<NumType>
          ([this, &col](const Trow & row) { return get_entry(row, col); });
    }
 
    DynList<Vector<Tcol, NumType>> to_rowlist() const
    {
      DynList<Vector<Tcol, NumType>> ret_val;
      DynList<Tcol> cols = col_domain_list();
      row_domain_list().
          for_each([this, &ret_val, cols](const Trow & row)
                     {
                       Vector<Tcol, NumType> vec(*cdomain_ptr, epsilon);
                       cols.for_each([this, &vec, &row](const Tcol & col)
                                       {
                                         auto *val_ptr =
                                             entries.search(RCPair(row, col));
                                         if (val_ptr != nullptr)
                                           vec.set_entry(col, val_ptr->second);
                                       });
                       ret_val.append(vec);
                     });
      return ret_val;
    }
 
    DynList<ColVector> to_collist() const
    {
      DynList<ColVector> ret_val;
      DynList<Trow> rows = row_domain_list();
      col_domain_list().
          for_each([this, &ret_val, &rows](const Tcol & col)
                     {
                       ColVector vec(*rdomain_ptr, epsilon);
                       rows.for_each([this, &vec, &col](const Trow & row)
                                       {
                                         auto *ptr = entries.search(RCPair(row, col));
                                         if (ptr != nullptr)
                                           vec.set_entry(row, ptr->second);
                                       });
                       ret_val.append(vec);
                     });
      return ret_val;
    }
 
    Vector<Tcol, NumType> get_row_vector(const Trow & row) const
    {
      ah_domain_error_if(not rdomain_ptr->has(row))
        << "get_row_vector(): invalid row";
      Vector<Tcol, NumType> ret_val(*cdomain_ptr);
      cdomain_ptr->for_each([&](const Tcol & col)
                              {
                                auto *ptr = entries.search(RCPair(row, col));
                                if (ptr != nullptr)
                                  ret_val.set_entry(col, ptr->second);
                              });
      return ret_val;
    }
 
    Vector<Trow, NumType> get_col_vector(const Tcol & col) const
    {
      ah_domain_error_if(not cdomain_ptr->has(col))
        << "get_col_vector(): invalid col";
      Vector<Trow, NumType> ret_val(*rdomain_ptr);
      rdomain_ptr->for_each([&](const Trow & row)
                              {
                                auto *ptr = entries.search(RCPair(row, col));
                                if (ptr != nullptr)
                                  ret_val.set_entry(row, ptr->second);
                              });
      return ret_val;
    }
 
    Vector<Trow, NumType>
    mult_matrix_vector_linear_comb(const Vector<Tcol, NumType> & vec) const
    {
      ah_domain_error_if(&vec.get_domain() != &get_col_domain())
        << "mult_matrix_vector_linear_comb: domain mismatch";
      Vector<Trow, NumType> ret_val(*rdomain_ptr);
      vec.for_each([this, &ret_val](const std::pair<Tcol, NumType> & p)
                     {
                       const Tcol & col = p.first;
                       const NumType & val = p.second;
                       ret_val += val * get_col_vector(col);
                     });
      return ret_val;
    }
 
    Vector<Trow, NumType>
    mult_matrix_vector_sparse(const Vector<Tcol, NumType> & vec) const
    {
      ah_domain_error_if(&vec.get_domain() != &get_col_domain())
        << "mult_matrix_vector_sparse: domain mismatch";
      Vector<Trow, NumType> ret_val(*rdomain_ptr);
      entries.for_each([&ret_val, &vec](const Pair & p)
                         {
                           const Trow & row = p.first.first;
                           const Tcol & col = p.first.second;
                           const NumType & entry_val = p.second;
                           NumType *vec_entry = vec.search_entry(col);
                           if (vec_entry != nullptr)
                             ret_val.set_entry(row, ret_val.get_entry(row) +
                                                    entry_val * *vec_entry);
                         });
      return ret_val;
    }
 
    Vector<Tcol, NumType>
    mult_vector_matrix_linear_comb(const Vector<Trow, NumType> & vec) const
    {
      ah_domain_error_if(&vec.get_domain() != &get_row_domain())
        << "mult_vector_matrix_linear_comb: domain mismatch";
      Vector<Tcol, NumType> ret_val(*cdomain_ptr);
      vec.for_each([this, &ret_val](const std::pair<Trow, NumType> & p)
                     {
                       const Trow & row = p.first;
                       const NumType & val = p.second;
                       ret_val += val * get_row_vector(row);
                     });
      return ret_val;
    }
 
    Vector<Trow, NumType>
    mult_matrix_vector_dot_product(const Vector<Tcol, NumType> & vec) const
    {
      ah_domain_error_if(&vec.get_domain() != &get_col_domain())
        << "mult_matrix_vector_dot_product: domain mismatch";
      Vector<Trow, NumType> ret_val(*rdomain_ptr);
      rdomain_ptr->for_each([&](const Trow & row)
                              {
                                ret_val.set_entry(row, vec * get_row_vector(row));
                              });
      return ret_val;
    }
 
    Vector<Tcol, NumType>
    mult_vector_matrix_dot_product(const Vector<Trow, NumType> & vec) const
    {
      ah_domain_error_if(&vec.get_domain() != &get_row_domain())
        << "mult_vector_matrix_dot_product: domain mismatch";
      Vector<Tcol, NumType> ret_val(*cdomain_ptr);
      cdomain_ptr->for_each([&](const Tcol & col)
                              {
                                ret_val.set_entry(col, vec * get_col_vector(col));
                              });
      return ret_val;
    }
 
    Vector<Trow, NumType>
    operator *(const Vector<Tcol, NumType> & vec) const
    {
      return mult_matrix_vector_linear_comb(vec);
    }
 
  private:
    template <typename T2row, typename T2col>
    void validate_domains_for_mult(const Matrix<T2row, T2col, NumType> & m)
    const
    {
      ah_domain_error_if(cdomain_ptr != m.rdomain_ptr)
        << "AxB: Column domain of A must be identical to row domain of B";
    }
 
  public:
    Matrix &set_vector_as_row(const Trow & row, const Vector<Tcol, NumType> & vec)
    {
      ah_domain_error_if(&vec.get_domain() != cdomain_ptr.get())
        << "set_vector_as_row: domain vec mismatch";
      vec.for_each([this, &row](const std::pair<Tcol, NumType> & p)
                     {
                       const Tcol & col = p.first;
                       const NumType & val = p.second;
                       set_entry(row, col, val);
                     });
      return *this;
    }
 
    Matrix &set_vector_as_col(const Tcol & col, const Vector<Trow, NumType> & vec)
    {
      ah_domain_error_if(&vec.get_domain() != rdomain_ptr.get())
        << "set_vector_as_col: domain vec mismatch";
      vec.for_each([this, &col](const std::pair<Trow, NumType> & p)
                     {
                       const Trow & row = p.first;
                       const NumType & val = p.second;
                       set_entry(row, col, val);
                     });
      return *this;
    }
 
    Matrix identity() const
    {
      static_assert(std::is_same_v<Trow, Tcol>,
                    "Row domain is not equal to column domain");
 
      Matrix<Trow, Tcol, NumType> i(*rdomain_ptr, *rdomain_ptr);
      rdomain_ptr->for_each([&i](const Trow & row)
                              {
                                i.set_entry(row, row, 1);
                              });
      return i;
    }
 
    Matrix &operator +=(const Matrix & m)
    {
      ah_domain_error_if(rdomain_ptr != m.rdomain_ptr or cdomain_ptr != m.cdomain_ptr)
        << "Matrix addition domain mismatch";
      m.entries.for_each([this](const Pair & p)
                           {
                             const Trow & row = p.first.first;
                             const Tcol & col = p.first.second;
                             const RCPair rcp = RCPair(row, col);
                             auto *ptr = entries.search(rcp);
                             if (ptr == nullptr)
                               entries.insert(rcp, p.second);
                             else
                               ptr->second += p.second;
                           });
      return *this;
    }
 
    Matrix operator +(const Matrix & m) const
    {
      Matrix ret_val = *this;
      ret_val += m;
      return ret_val;
    }
 
    Matrix &operator -=(const Matrix & m)
    {
      ah_domain_error_if(rdomain_ptr != m.rdomain_ptr or cdomain_ptr != m.cdomain_ptr)
        << "Matrix subtraction domain mismatch";
      m.entries.for_each([this](const Pair & p)
                           {
                             const Trow & row = p.first.first;
                             const Tcol & col = p.first.second;
                             const RCPair rcp = RCPair(row, col);
                             auto *ptr = entries.search(rcp);
                             if (ptr == nullptr)
                               entries.insert(rcp, -p.second);
                             else
                               ptr->second -= p.second;
                           });
      return *this;
    }
 
    Matrix operator -(const Matrix & m) const
    {
      Matrix ret_val = *this;
      ret_val -= m;
      return ret_val;
    }
 
    Matrix &mult_by_scalar(const NumType & scalar)
    {
      entries.for_each([&scalar](const Pair & p)
                         {
                           NumType & val = const_cast<Pair &>(p).second;
                           val *= scalar;
                         });
      return *this;
    }
 
    template <typename T2row, typename T2col>
    Matrix<Trow, T2col>
    vector_matrix_mult(const Matrix<T2row, T2col, NumType> & m) const
    {
      validate_domains_for_mult(m);
      Matrix<Trow, T2col> ret_val(*rdomain_ptr, m.get_col_domain());
      rdomain_ptr->
          for_each([this, &m, &ret_val](const Trow & row)
                     {
                       ret_val.set_vector_as_row(row, get_row_vector(row) * m);
                     });
      return ret_val;
    }
 
    template <typename T2row, typename T2col>
    Matrix<Trow, T2col>
    matrix_vector_mult(const Matrix<T2row, T2col, NumType> & m) const
    {
      validate_domains_for_mult(m);
      Matrix<Trow, T2col> ret_val(*rdomain_ptr, m.get_col_domain());
      m.get_col_domain().
        for_each([this, &m, &ret_val](const Tcol & col)
                   {
                     ret_val.set_vector_as_col(col,
                                               *this * m.get_col_vector(col));
                   });
      return ret_val;
    }
 
    void print() const
    {
      std::cout << "Non zero entries = " << '\n';
      entries.for_each([](const Pair & p)
                         {
                           std::cout << "(" << Aleph::to_str(p.first.first) << ","
                               << Aleph::to_str(p.first.second) << ") : "
                               << Aleph::to_str(p.second) << '\n';
                         });
      std::cout << '\n';
    }
 
    [[nodiscard]] std::string to_str() const
    {
      typedef std::pair<std::string, std::string> StrPair;
      DynList<ColVector> cols = to_collist();
      // list of longest strings for each column
      const DynList<std::string> collens = cols.
          template maps<std::string>([](const ColVector & vec)
                                       {
                                         return vec.to_list().
                                                    template foldl<std::string>
                                                    ("", [](const std::string & s, const NumType & val)
                                                       {
                                                         std::string val_str = Aleph::to_str(val);
                                                         return val_str.size() > s.size() ? val_str : s;
                                                       });
                                       });
 
      const DynList<std::string> col_labels =
          sort(cdomain_ptr->template maps<std::string>([](const Tcol & col)
                                                         {
                                                           return Aleph::to_str(col);
                                                         }));
 
      DynList<std::string> longers = zip(collens, col_labels).
          template maps<std::string>([](const StrPair & p)
                                       {
                                         return p.first.size() > p.second.size() ?
                                                  p.first :
                                                  p.second;
                                       });
 
      DynList<std::string> row_labels =
          sort(rdomain_ptr->template maps<std::string>([](const Trow & row)
                                                         {
                                                           return Aleph::to_str(row);
                                                         }));
 
      const std::string longer_row_label = row_labels.template foldl<std::string>
          (std::string(), [](const std::string & s,
                             const std::string & row_label)
             {
               return row_label.size() > s.size() ? row_label : s;
             });
 
      size_t longer_row_label_sz = longer_row_label.size();
 
      std::string ret_val = std::string(longer_row_label_sz, ' ') + " |";
      ret_val = zip(col_labels, longers).template foldl<std::string>
          (ret_val, [](const std::string & s, const StrPair & p)
             {
               const size_t fsz = p.first.size();
               const size_t ssz = p.second.size();
               return s + std::string((fsz > ssz ? fsz - ssz : ssz - fsz) + 1, ' ') + p.first;
             });
 
      ret_val += "\n" + std::string(ret_val.size(), '-');
 
      zip(row_labels, row_domain_list()).for_each
          ([this, &ret_val, &longers, &longer_row_label_sz]
       (const std::pair<std::string, Trow> & p)
             {
               const std::string & label = p.first;
               const Trow & row = p.second;
 
               ret_val += "\n" +
                   std::string(longer_row_label_sz - label.size(), ' ') + label + " |";
               ret_val =
                   zip(get_row_as_list(row), longers).template foldl<std::string>
                   (ret_val, [](const std::string & s,
                                const std::pair<NumType, std::string> & q)
                      {
                        const std::string val = Aleph::to_str(q.first);
                        const std::string & len_str = q.second;
                        assert(len_str.size() - val.size() + 1);
                        return s + std::string(len_str.size() - val.size() + 1, ' ') + val;
                      });
             });
 
      return ret_val;
    }
 
    bool are_equal(const NumType & n1, const NumType & n2) const
    {
      return is_zero(n1 - n2);
    }
 
    bool equal_to(const Matrix & m) const
    {
      ah_domain_error_if(not (rdomain_ptr == m.rdomain_ptr and cdomain_ptr == m.cdomain_ptr))
        << "invalid matrix's domains";
      return entries.all([&m, this](const Pair & p)
                           {
                             const Trow & row = p.first.first;
                             const Tcol & col = p.first.second;
                             return are_equal(m.get_entry(row, col), p.second);
                           }) and
             m.entries.all([this](const Pair & p)
                             {
                               const Trow & row = p.first.first;
                               const Tcol & col = p.first.second;
                               return are_equal(get_entry(row, col), p.second);
                             });
    }
 
    bool operator ==(const Matrix & m) const { return equal_to(m); }
 
    bool operator !=(const Matrix & m) const { return not equal_to(m); }
  };
 
  template <typename Trow, typename Tcol, typename NumType>
  inline
  Matrix<Trow, Tcol, NumType> operator *(const NumType & scalar,
                                         const Matrix<Trow, Tcol, NumType> & m)
  {
    Matrix<Trow, Tcol, NumType> ret_val = m;
    return ret_val.mult_by_scalar(scalar);
  }
 
  template <typename Trow, typename Tcol, typename NumType>
  inline
  Vector<Tcol, NumType> operator *(const Vector<Trow, NumType> & vec,
                                   const Matrix<Trow, Tcol, NumType> & m)
  {
    return m.mult_vector_matrix_linear_comb(vec);
  }
 
  template <typename Trow, typename Tcol, typename NumType>
  const NumType Matrix<Trow, Tcol, NumType>::default_epsilon = 1e-7;
 
  template <typename Trow, typename Tcol, typename NumType>
  inline
  std::ostream &operator <<(std::ostream & s, const Matrix<Trow, Tcol, NumType> & m)
  {
    return s << m.to_str();
  }
 
  template <typename Trow, typename Tcol, typename NumType = double>
  inline
  Matrix<Trow, Tcol, NumType> outer_product(const Vector<Trow, NumType> & v1,
                                            const Vector<Tcol, NumType> & v2)
  {
    Matrix<Trow, Tcol, NumType> ret_val(v1.get_domain(), v2.get_domain());
 
    v1.for_each([&ret_val, &v2](const std::pair<Trow, NumType> & p1)
                  {
                    v2.for_each([&ret_val, &p1](const std::pair<Tcol, NumType> & p2)
                                  {
                                    ret_val.set_entry(p1.first, p2.first,
                                                      p1.second * p2.second);
                                  });
                  });
    return ret_val;
  }
} // end namespace Aleph
 
# endif // AH_MATRIX_H