densematrix.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
/*
 * File:   densematrix.hh
 * Author: ngo
 *
 * Created on April 15, 2011
 */
 
#ifndef DENSEMATRIX_HH
#define DENSEMATRIX_HH
 
#include <cstdlib>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <string>
 
#include "exceptions.hh"
#include "sparsematrix.hh"
#include "vector.hh"
 
namespace hdnum {
 
// forward-declare the sparse matrix template to make the transforming
// constructor from hdnum::SparseMatrix -> hdnum::DenseMatrix working
template <typename REAL>
class SparseMatrix;
 
template <typename REAL>
class DenseMatrix {
public:
    typedef std::size_t size_type;
    typedef typename std::vector<REAL> VType;
    typedef typename VType::const_iterator ConstVectorIterator;
    typedef typename VType::iterator VectorIterator;
 
private:
    VType m_data;        // Matrix data is stored in an STL vector!
    std::size_t m_rows;  // Number of Matrix rows
    std::size_t m_cols;  // Number of Matrix columns
 
    static bool bScientific;
    static std::size_t nIndexWidth;
    static std::size_t nValueWidth;
    static std::size_t nValuePrecision;
 
    REAL myabs(REAL x) const {
        if (x >= REAL(0))
            return x;
        else
            return -x;
    }
 
    inline REAL& at(const std::size_t row, const std::size_t col) {
        return m_data[row * m_cols + col];
    }
 
    inline const REAL& at(const std::size_t row, const std::size_t col) const {
        return m_data[row * m_cols + col];
    }
 
public:
    DenseMatrix() : m_data(0, 0), m_rows(0), m_cols(0) {}
 
    DenseMatrix(const std::size_t _rows, const std::size_t _cols,
                const REAL def_val = 0)
        : m_data(_rows * _cols, def_val), m_rows(_rows), m_cols(_cols) {}
 
    DenseMatrix(const std::initializer_list<std::initializer_list<REAL>>& v) {
        m_rows = v.size();
        m_cols = v.begin()->size();
        for (auto row : v) {
            if (row.size() != m_cols) {
                std::cout << "Zeilen der Matrix nicht gleich lang" << std::endl;
                exit(1);
            }
            for (auto elem : row) m_data.push_back(elem);
        }
    }
 
    DenseMatrix(const hdnum::SparseMatrix<REAL>& other)
        : m_data(other.rowsize() * other.colsize()), m_rows(other.rowsize()),
          m_cols(other.colsize()) {
        using counter_type = typename hdnum::SparseMatrix<REAL>::size_type;
        counter_type row_index {};
        for (auto& row : other) {
            for (auto it = row.ibegin(); it != row.iend(); it++) {
                this->operator[](row_index)[it.index()] = it.value();
            }
            row_index++;
        }
    }
 
    void addNewRow(const hdnum::Vector<REAL>& rowvector) {
        m_rows++;
        m_cols = rowvector.size();
        for (std::size_t i = 0; i < m_cols; i++) m_data.push_back(rowvector[i]);
    }
 
    /*
    // copy constructor (not needed, since it inherits from the STL vector)
    DenseMatrix( const DenseMatrix& A )
    {
    this->m_data = A.m_data;
    m_rows = A.m_rows;
    m_cols = A.m_cols;
    }
    */
 
    size_t rowsize() const { return m_rows; }
 
    size_t colsize() const { return m_cols; }
 
    // pretty-print output properties
    bool scientific() const { return bScientific; }
 
    void scientific(bool b) const { bScientific = b; }
 
    std::size_t iwidth() const { return nIndexWidth; }
 
    std::size_t width() const { return nValueWidth; }
 
    std::size_t precision() const { return nValuePrecision; }
 
    void iwidth(std::size_t i) const { nIndexWidth = i; }
 
    void width(std::size_t i) const { nValueWidth = i; }
 
    void precision(std::size_t i) const { nValuePrecision = i; }
 
    // overloaded element access operators
    // write access on matrix element A_ij using A(i,j)
    inline REAL& operator()(const std::size_t row, const std::size_t col) {
        assert(row < m_rows || col < m_cols);
        return at(row, col);
    }
 
    inline const REAL& operator()(const std::size_t row,
                                  const std::size_t col) const {
        assert(row < m_rows || col < m_cols);
        return at(row, col);
    }
 
    const ConstVectorIterator operator[](const std::size_t row) const {
        assert(row < m_rows);
        return m_data.begin() + row * m_cols;
    }
 
    VectorIterator operator[](const std::size_t row) {
        assert(row < m_rows);
        return m_data.begin() + row * m_cols;
    }
 
    DenseMatrix& operator=(const DenseMatrix& A) {
        m_data = A.m_data;
        m_rows = A.m_rows;
        m_cols = A.m_cols;
        return *this;
    }
 
    DenseMatrix& operator=(const REAL value) {
        for (std::size_t i = 0; i < rowsize(); i++)
            for (std::size_t j = 0; j < colsize(); j++) (*this)(i, j) = value;
        return *this;
    }
 
    DenseMatrix sub(size_type i, size_type j, size_type rows, size_type cols) {
        DenseMatrix A(rows, cols);
        DenseMatrix& self = *this;
        for (size_type k1 = 0; k1 < rows; k1++) {
            for (size_type k2 = 0; k2 < cols; k2++) {
                A[k1][k2] = self[k1 + i][k2 + j];
            }
        }
        return A;
    }
 
    DenseMatrix transpose() const {
        DenseMatrix A(m_cols, m_rows);
        for (size_type i = 0; i < m_rows; i++) {
            for (size_type j = 0; j < m_cols; j++) {
                A[j][i] = this->operator[](i)[j];
            }
        }
        return A;
    }
 
    // Basic Matrix Operations
 
    DenseMatrix& operator+=(const DenseMatrix& B) {
        for (size_type i = 0; i < rowsize(); ++i) {
            for (size_type j = 0; j < colsize(); ++j) {
                (*this)(i, j) += B(i, j);
            }
        }
        return *this;
    }
 
    DenseMatrix& operator-=(const DenseMatrix& B) {
        for (std::size_t i = 0; i < rowsize(); ++i)
            for (std::size_t j = 0; j < colsize(); ++j)
                (*this)(i, j) -= B(i, j);
        return *this;
    }
 
    DenseMatrix& operator*=(const REAL s) {
        for (std::size_t i = 0; i < rowsize(); ++i)
            for (std::size_t j = 0; j < colsize(); ++j) (*this)(i, j) *= s;
        return *this;
    }
 
    DenseMatrix& operator/=(const REAL s) {
        for (std::size_t i = 0; i < rowsize(); ++i)
            for (std::size_t j = 0; j < colsize(); ++j) (*this)(i, j) /= s;
        return *this;
    }
 
    void update(const REAL s, const DenseMatrix& B) {
        for (std::size_t i = 0; i < rowsize(); ++i)
            for (std::size_t j = 0; j < colsize(); ++j)
                (*this)(i, j) += s * B(i, j);
    }
 
    template <class V>
    void mv(Vector<V>& y, const Vector<V>& x) const {
        if (this->rowsize() != y.size())
            HDNUM_ERROR("mv: size of A and y do not match");
        if (this->colsize() != x.size())
            HDNUM_ERROR("mv: size of A and x do not match");
        for (std::size_t i = 0; i < rowsize(); ++i) {
            y[i] = 0;
            for (std::size_t j = 0; j < colsize(); ++j)
                y[i] += (*this)(i, j) * x[j];
        }
    }
 
    template <class V>
    void umv(Vector<V>& y, const Vector<V>& x) const {
        if (this->rowsize() != y.size())
            HDNUM_ERROR("mv: size of A and y do not match");
        if (this->colsize() != x.size())
            HDNUM_ERROR("mv: size of A and x do not match");
        for (std::size_t i = 0; i < rowsize(); ++i) {
            for (std::size_t j = 0; j < colsize(); ++j)
                y[i] += (*this)(i, j) * x[j];
        }
    }
 
    template <class V>
    void umv(Vector<V>& y, const V& s, const Vector<V>& x) const {
        if (this->rowsize() != y.size())
            HDNUM_ERROR("mv: size of A and y do not match");
        if (this->colsize() != x.size())
            HDNUM_ERROR("mv: size of A and x do not match");
        for (std::size_t i = 0; i < rowsize(); ++i) {
            for (std::size_t j = 0; j < colsize(); ++j)
                y[i] += s * (*this)(i, j) * x[j];
        }
    }
 
    void mm(const DenseMatrix<REAL>& A, const DenseMatrix<REAL>& B) {
        if (this->rowsize() != A.rowsize())
            HDNUM_ERROR("mm: size incompatible");
        if (this->colsize() != B.colsize())
            HDNUM_ERROR("mm: size incompatible");
        if (A.colsize() != B.rowsize()) HDNUM_ERROR("mm: size incompatible");
 
        for (std::size_t i = 0; i < rowsize(); i++)
            for (std::size_t j = 0; j < colsize(); j++) {
                (*this)(i, j) = 0;
                for (std::size_t k = 0; k < A.colsize(); k++)
                    (*this)(i, j) += A(i, k) * B(k, j);
            }
    }
 
    void umm(const DenseMatrix<REAL>& A, const DenseMatrix<REAL>& B) {
        if (this->rowsize() != A.rowsize())
            HDNUM_ERROR("mm: size incompatible");
        if (this->colsize() != B.colsize())
            HDNUM_ERROR("mm: size incompatible");
        if (A.colsize() != B.rowsize()) HDNUM_ERROR("mm: size incompatible");
 
        for (std::size_t i = 0; i < rowsize(); i++)
            for (std::size_t j = 0; j < colsize(); j++)
                for (std::size_t k = 0; k < A.colsize(); k++)
                    (*this)(i, j) += A(i, k) * B(k, j);
    }
 
    void sc(const Vector<REAL>& x, std::size_t k) {
        if (this->rowsize() != x.size()) HDNUM_ERROR("cc: size incompatible");
 
        for (std::size_t i = 0; i < rowsize(); i++) (*this)(i, k) = x[i];
    }
 
    void sr(const Vector<REAL>& x, std::size_t k) {
        if (this->colsize() != x.size()) HDNUM_ERROR("cc: size incompatible");
 
        for (std::size_t i = 0; i < colsize(); i++) (*this)(k, i) = x[i];
    }
 
    REAL norm_infty() const {
        REAL norm(0.0);
        for (std::size_t i = 0; i < rowsize(); i++) {
            REAL sum(0.0);
            for (std::size_t j = 0; j < colsize(); j++)
                sum += myabs((*this)(i, j));
            if (sum > norm) norm = sum;
        }
        return norm;
    }
 
    REAL norm_1() const {
        REAL norm(0.0);
        for (std::size_t j = 0; j < colsize(); j++) {
            REAL sum(0.0);
            for (std::size_t i = 0; i < rowsize(); i++)
                sum += myabs((*this)(i, j));
            if (sum > norm) norm = sum;
        }
        return norm;
    }
 
    Vector<REAL> operator*(const Vector<REAL>& x) const {
        assert(x.size() == colsize());
 
        Vector<REAL> y(rowsize());
        for (std::size_t r = 0; r < rowsize(); ++r) {
            for (std::size_t c = 0; c < colsize(); ++c) {
                y[r] += at(r, c) * x[c];
            }
        }
        return y;
    }
 
    DenseMatrix operator*(const DenseMatrix& x) const {
        assert(colsize() == x.rowsize());
 
        const std::size_t out_rows = rowsize();
        const std::size_t out_cols = x.colsize();
        DenseMatrix y(out_rows, out_cols, 0.0);
        for (std::size_t r = 0; r < out_rows; ++r)
            for (std::size_t c = 0; c < out_cols; ++c)
                for (std::size_t i = 0; i < colsize(); ++i)
                    y(r, c) += at(r, i) * x(i, c);
 
        return y;
    }
 
    DenseMatrix operator+(const DenseMatrix& x) const {
        assert(colsize() == x.colsize());
        assert(rowsize() == x.rowsize());
 
        const std::size_t out_rows = rowsize();
        const std::size_t out_cols = x.colsize();
        DenseMatrix y(out_rows, out_cols, 0.0);
        y = *this;
        y += x;
        return y;
    }
 
    DenseMatrix operator-(const DenseMatrix& x) const {
        assert(colsize() == x.colsize());
        assert(rowsize() == x.rowsize());
 
        const std::size_t out_rows = rowsize();
        const std::size_t out_cols = x.colsize();
        DenseMatrix y(out_rows, out_cols, 0.0);
        y = *this;
        y -= x;
        return y;
    }
};
 
template <typename REAL>
bool DenseMatrix<REAL>::bScientific = true;
template <typename REAL>
std::size_t DenseMatrix<REAL>::nIndexWidth = 10;
template <typename REAL>
std::size_t DenseMatrix<REAL>::nValueWidth = 10;
template <typename REAL>
std::size_t DenseMatrix<REAL>::nValuePrecision = 3;
 
template <typename REAL>
inline std::ostream& operator<<(std::ostream& s, const DenseMatrix<REAL>& A) {
    s << std::endl;
    s << " " << std::setw(A.iwidth()) << " "
      << "  ";
    for (typename DenseMatrix<REAL>::size_type j = 0; j < A.colsize(); ++j)
        s << std::setw(A.width()) << j << " ";
    s << std::endl;
 
    for (typename DenseMatrix<REAL>::size_type i = 0; i < A.rowsize(); ++i) {
        s << " " << std::setw(A.iwidth()) << i << "  ";
        for (typename DenseMatrix<REAL>::size_type j = 0; j < A.colsize();
             ++j) {
            if (A.scientific()) {
                s << std::setw(A.width()) << std::scientific << std::showpoint
                  << std::setprecision(A.precision()) << A[i][j] << " ";
            } else {
                s << std::setw(A.width()) << std::fixed << std::showpoint
                  << std::setprecision(A.precision()) << A[i][j] << " ";
            }
        }
        s << std::endl;
    }
    return s;
}
 
template <typename REAL>
inline void fill(DenseMatrix<REAL>& A, const REAL& t) {
    for (typename DenseMatrix<REAL>::size_type i = 0; i < A.rowsize(); ++i)
        for (typename DenseMatrix<REAL>::size_type j = 0; j < A.colsize(); ++j)
            A[i][j] = t;
}
 
template <typename REAL>
inline void zero(DenseMatrix<REAL>& A) {
    for (std::size_t i = 0; i < A.rowsize(); ++i)
        for (std::size_t j = 0; j < A.colsize(); ++j) A(i, j) = REAL(0);
}
 
template <class T>
inline void identity(DenseMatrix<T>& A) {
    for (typename DenseMatrix<T>::size_type i = 0; i < A.rowsize(); ++i)
        for (typename DenseMatrix<T>::size_type j = 0; j < A.colsize(); ++j)
            if (i == j)
                A[i][i] = T(1);
            else
                A[i][j] = T(0);
}
 
template <typename REAL>
inline void spd(DenseMatrix<REAL>& A) {
    if (A.rowsize() != A.colsize() || A.rowsize() == 0)
        HDNUM_ERROR("need square and nonempty matrix");
    for (std::size_t i = 0; i < A.rowsize(); ++i)
        for (std::size_t j = 0; j < A.colsize(); ++j)
            if (i == j)
                A(i, i) = REAL(4.0);
            else
                A(i, j) = -REAL(1.0) / ((i - j) * (i - j));
}
 
template <typename REAL>
inline void vandermonde(DenseMatrix<REAL>& A, const Vector<REAL> x) {
    if (A.rowsize() != A.colsize() || A.rowsize() == 0)
        HDNUM_ERROR("need square and nonempty matrix");
    if (A.rowsize() != x.size()) HDNUM_ERROR("need A and x of same size");
    for (typename DenseMatrix<REAL>::size_type i = 0; i < A.rowsize(); ++i) {
        REAL p(1.0);
        for (typename DenseMatrix<REAL>::size_type j = 0; j < A.colsize();
             ++j) {
            A[i][j] = p;
            p *= x[i];
        }
    }
}
 
template <typename REAL>
inline void gnuplot(const std::string& fname, const DenseMatrix<REAL>& A) {
    std::fstream f(fname.c_str(), std::ios::out);
    for (typename DenseMatrix<REAL>::size_type i = 0; i < A.rowsize(); ++i) {
        for (typename DenseMatrix<REAL>::size_type j = 0; j < A.colsize();
             ++j) {
            if (A.scientific()) {
                f << std::setw(A.width()) << std::scientific << std::showpoint
                  << std::setprecision(A.precision()) << A[i][j];
            } else {
                f << std::setw(A.width()) << std::fixed << std::showpoint
                  << std::setprecision(A.precision()) << A[i][j];
            }
        }
        f << std::endl;
    }
    f.close();
}
 
template <typename REAL>
inline void readMatrixFromFileDat(const std::string& filename,
                                  DenseMatrix<REAL>& A) {
    std::string buffer;
    std::ifstream fin(filename.c_str());
    std::size_t i = 0;
    std::size_t j = 0;
    if (fin.is_open()) {
        while (std::getline(fin, buffer)) {
            std::istringstream iss(buffer);
            hdnum::Vector<REAL> rowvector;
            while (iss) {
                std::string sub;
                iss >> sub;
                // std::cout << " sub = " << sub.c_str() << ": ";
                if (sub.length() > 0) {
                    REAL a = atof(sub.c_str());
                    // std::cout << std::fixed << std::setw(10) <<
                    // std::setprecision(5) << a;
                    rowvector.push_back(a);
                }
                j++;
            }
            if (rowvector.size() > 0) {
                A.addNewRow(rowvector);
                i++;
                // std::cout << std::endl;
            }
        }
        fin.close();
    } else {
        HDNUM_ERROR("Could not open file!");
    }
}
 
template <typename REAL>
inline void readMatrixFromFileMatrixMarket(const std::string& filename,
                                           DenseMatrix<REAL>& A) {
    std::string buffer;
    std::ifstream fin(filename.c_str());
    std::size_t i = 0;
    std::size_t j = 0;
    if (fin.is_open()) {
        // ignore all comments from the file (starting with %)
        while (fin.peek() == '%') fin.ignore(2048, '\n');
 
        std::getline(fin, buffer);
        std::istringstream first_line(buffer);
        first_line >> i >> j;
        DenseMatrix<REAL> A_temp(i, j);
 
        while (std::getline(fin, buffer)) {
            std::istringstream iss(buffer);
 
            REAL value {};
            iss >> i >> j >> value;
            // i-1, j-1, because matrix market does not use zero based indexing
            A_temp(i - 1, j - 1) = value;
        }
        A = A_temp;
        fin.close();
    } else {
        HDNUM_ERROR("Could not open file! \"" + filename + "\"");
    }
}
 
}  // namespace hdnum
 
#endif  // DENSEMATRIX_HH