matrix.h

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 *                                                                           *
 *                       / /   / /    __    \ \  / /                         *
 *                      / /__ / /   / _ \    \ \/ /                          *
 *                     /  ___  /   | |/_/    / /\ \                          *
 *                    / /   / /    \_\      / /  \ \                         *
 *                                                                           *
 *                         Jakub Benda (c) 2014                              *
 *                     Charles University in Prague                          *
 *                                                                           *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#ifndef HEX_SPMATRIX
#define HEX_SPMATRIX

#include <cassert>
#include <complex>
#include <cstring>
#include <iostream>
#include <vector>

#ifndef NO_PNG
#include <png++/png.hpp>
#endif

#include <umfpack.h>

#include "arrays.h"
#include "complex.h"

// forward declaration of classes in order to enable them as return types
// of other classes before proper definition
template <class T> class RowMatrix;
template <class T> class ColMatrix;
class CooMatrix;
class CscMatrix;
class CsrMatrix;
class SymDiaMatrix;

template <class T> class DenseMatrix
{
    public:
        
        // constructors
        DenseMatrix()
            : rows_(0), cols_(0), data_() {}
        DenseMatrix(int rows, int cols)
            : rows_(rows), cols_(cols), data_(rows * cols) {}
        DenseMatrix(int rows, int cols, const ArrayView<T> data)
            : rows_(rows), cols_(cols), data_(data) { assert(data.size() == (size_t)(rows * cols)); }
        
        // explicit conversion to cArrayView
        ArrayView<T> data () { return data_; }
        const ArrayView<T> data () const { return data_; }
        
        // getters
        size_t size () const { return rows_ * cols_; }
        int cols () const { return cols_; }
        int rows () const { return rows_; }
        T * begin () { return data_.begin(); }
        T const * begin () const { return data_.begin(); }
        
    private:
        
        int rows_;
        
        int cols_;
        
        NumberArray<T> data_;
};

template <class Type> class ColMatrix : public DenseMatrix<Type>
{
    public:
        
        // constructor
        ColMatrix ()
            : DenseMatrix<Type>() {}
        ColMatrix (int size)
            : DenseMatrix<Type>(size, size) {}
        ColMatrix (int rows, int cols)
            : DenseMatrix<Type>(rows, cols) {}
        ColMatrix (int rows, int cols, const ArrayView<Type> data)
            : DenseMatrix<Type>(rows, cols, data) {}
        
        explicit ColMatrix (DenseMatrix<Type> const & M)
            : DenseMatrix<Type>(M) {}
        
        template <class Function> void populate (Function f)
        {
            // pointer to data
            Type * ptr = this->begin();
            
            // fill the matrix
            for (int icol = 0; icol < this->cols(); icol++)
            for (int irow = 0; irow < this->rows(); irow++)
                *(ptr++) = f(irow,icol);
        }
        
        RowMatrix<Type> T () const
        {
            return RowMatrix<Type>
            (
                this->cols(),
                this->rows(),
                this->data()
            );
        }
        
        ArrayView<Type> col (int i)
        {
            return ArrayView<Type>
            (
                this->data(),
                i * this->rows(),
                this->rows()
            );
        }
        const ArrayView<Type> col (int i) const
        {
            return ArrayView<Type>
            (
                this->data(),
                i * this->rows(),
                this->rows()
            );
        }
        

        Type operator() (int i, int j) const { return col(j)[i]; }
        Type & operator() (int i, int j) { return col(j)[i]; }
};

template <class Type> class RowMatrix : public DenseMatrix<Type>
{
    public:
        
        // constructor
        RowMatrix ()
            : DenseMatrix<Type>() {}
        RowMatrix (int size)
            : DenseMatrix<Type>(size, size) {}
        RowMatrix (int rows, int cols)
            : DenseMatrix<Type>(rows, cols) {}
        RowMatrix (int rows, int cols, const ArrayView<Type> data)
            : DenseMatrix<Type>(rows, cols, data) {}
        RowMatrix (ColMatrix<Type> const & m)
            : DenseMatrix<Type>(m.rows(), m.cols(), m.data()) { reorder_(); }
        
        explicit RowMatrix (DenseMatrix<Type> const & M)
            : DenseMatrix<Type>(M) {}
        
        template <class Function> void populate (Function f)
        {
            // pointer to data
            Type * ptr = this->begin();
            
            // fill the matrix
            for (int irow = 0; irow < this->rows(); irow++)
            for (int icol = 0; icol < this->cols(); icol++)
                *(ptr++) = f(irow,icol);
        }
        
        template <class Function> void transform (Function f)
        {
            // pointer to data
            Type * ptr = this->begin();
            
            // transform the matrix
            for (int irow = 0; irow < this->rows(); irow++)
            for (int icol = 0; icol < this->cols(); icol++)
                f(irow,icol,*(ptr++));
        }
        
        ColMatrix<Type> T() const
        {
            return ColMatrix<Type>
            (
                this->cols(),
                this->rows(),
                this->data()
            );
        }
        
        ArrayView<Type> row (int i)
        {
            return ArrayView<Type>
            (
                this->data(),
                i * this->cols(),
                this->cols()
            );
        }
        const ArrayView<Type> row (int i) const
        {
            return ArrayView<Type>
            (
                this->data(),
                i * this->cols(),
                this->cols()
            );
        }
        

        Type operator() (int i, int j) const { return row(i)[j]; }
        Type & operator() (int i, int j) { return row(i)[j]; }
        
        Type operator() (const ArrayView<Type> u, const ArrayView<Type> v) const
        {
            // return value
            Type sum = 0;
            
            // pointer to matrix elements
            Type const * restrict pA = this->data().data();
            
            // pointers to vector elements
            Type const * const restrict pu = u.data();
            Type const * const restrict pv = v.data();
            
            // for all elements
            for (int i = 0; i < this->rows(); i++)
            for (int j = 0; j < this->cols(); j++)
                sum += pu[i] * (*(pA++)) * pv[j];
            
            return sum;
        }
        
        static RowMatrix<Type> Eye (int size)
        {
            RowMatrix<Type> M (size, size);
            
            for (int i = 0; i < 0; i++)
                M.data()[i * size + i] = Type(1);
            
            return M;
        }
        
        // arithmetic operators
        RowMatrix<Type> const & operator += (DenseMatrix<Type> const & A)
        {
            assert (this->rows() == A.rows());
            assert (this->cols() == A.cols());
            
            this->data() += A.data();
            
            return *this;
        }
        RowMatrix<Type> const & operator -= (DenseMatrix<Type> const & A)
        {
            assert (this->rows() == A.rows());
            assert (this->cols() == A.cols());
            
            this->data() -= A.data();
            
            return *this;
        }
        RowMatrix<Type> const & operator *= (Type x)
        {
            this->data() *= x;
            return *this;
        }
        
        void write (std::ostream & out, std::string const & pre = "", std::string const & pos = "") const
        {
            // data pointer
            Type const * ptr = this->data().begin();
            
            for (int irow = 0; irow < this->rows(); irow++)
            {
                out << pre;
                for (int icol = 0; icol < this->cols(); icol++)
                {
                    if (*ptr >= 0) out << " "; out << *ptr++ << " ";
                }
                
                out << pos << "\n";
            }
        }
        
#ifndef NO_PNG

        void plot_abs (std::ofstream & out) const
        {
            // create empty gray-scale 1-bit image
            png::image<png::gray_pixel_16> image (this->cols(), this->rows());
            
            // skip empty matrix
            if (this->data().size() > 0)
            {
                // find minimal and maximal value
                double min = std::abs(this->data()[0]), max = std::abs(this->data()[0]);
                for (Type const & x : this->data())
                {
                    double absx = std::abs(x);
                    if (absx < min) min = absx;
                    if (absx > max) max = absx;
                }
                
                // for all elements
                for (int y = 0; y < this->rows(); y++)
                {
                    for (int x = 0; x < this->cols(); x++)
                    {
                        double fraction = 1. - (std::abs(this->data()[y * this->cols() + x]) - min) / (max - min);
                        image.set_pixel(x, y, std::ceil(65535 * fraction));
                    }
                }
            }
            
            // save image
            image.write_stream(out);
        }
#endif
        
    private:
        
        void reorder_ ()
        {
            NumberArray<Type> new_data (this->data().size());
            
            for (int irow = 0; irow < this->rows(); irow++)
            for (int icol = 0; icol < this->cols(); icol++)
                new_data[irow * this->cols() + icol] = this->data()[icol * this->rows() + irow];
            
            this->data() = new_data;
        }
};

template <class T> RowMatrix<T> operator + (RowMatrix<T> const & A, RowMatrix<T> const & B) { RowMatrix<T> C(A); C.data() += B.data(); return C; }
template <class T> RowMatrix<T> operator - (RowMatrix<T> const & A, RowMatrix<T> const & B) { RowMatrix<T> C(A); C.data() -= B.data(); return C; }
template <class T> RowMatrix<T> operator * (T x, RowMatrix<T> const & A) { RowMatrix<T> B(A); B.data() *= x; return B; }
template <class T> RowMatrix<T> operator * (RowMatrix<T> const & A, T x) { RowMatrix<T> B(A); B.data() *= x; return B; }
template <class T> RowMatrix<T> operator / (RowMatrix<T> const & A, T x) { RowMatrix<T> B(A); B.data() /= x; return B; }

// matrix multiplication
template <class T> RowMatrix<T> operator * (RowMatrix<T> const & A, ColMatrix<T> const & B);

template <class Type> RowMatrix<Type> operator * (RowMatrix<Type> const & A, ColMatrix<Type> const & B)
{
    assert(A.cols() == B.rows());
    
    // sizes
    int rows = A.rows();
    int cols = B.cols();
    int comm = A.cols();
    
    // output matrix, initialized to zero
    RowMatrix<Type> C(A.rows(), B.cols());
    
    // data pointers
    Type const * restrict pA = A.begin();
    Type const * restrict pB = B.begin();
    Type       * restrict pC = C.begin();
    
    // for all rows of A
    for (int irow = 0; irow < rows; irow++)
    {
        // for all columns of B
        for (int icol = 0; icol < cols; icol++)
        {
            // compute the scalar product "row * column"
            for (int k = 0; k < comm; k++)
                (*pC) += (*(pA + k)) * (*pB++);
            
            // move to next element of C
            pC++;
        }
        
        // move to the next row of A
        pA += A.cols();
        
        // reset B's data pointer
        pB = B.begin();
    }
    
    // return result
    return C;
}

typedef enum {
    none         = 0,  // b000,
    strict_lower = 1,  // b001, strict_lower
    strict_upper = 2,  // b010,                           strict_upper
    strict_both  = 3,  // b011, strict_lower |            strict_upper
    diagonal     = 4,  // b100,                diagonal
    lower        = 5,  // b101, strict_lower | diagonal
    upper        = 6,  // b110,                diagonal | strict_upper
    both         = 7   // b111, strict_lower | diagonal | strict_upper
} MatrixTriangle;

class CscMatrix
{
    
public:
    
    // Constructors
    
    CscMatrix()
        : m_(0), n_(0) {}
    CscMatrix(size_t m, size_t n)
        : m_(m), n_(n) {}
    CscMatrix(CscMatrix const & A)
        : m_(A.m_), n_(A.n_), p_(A.p_), i_(A.i_), x_(A.x_) {}
    CscMatrix(size_t m, size_t n, const lArrayView p, const lArrayView i, const cArrayView x)
        : m_(m), n_(n), p_(p), i_(i), x_(x) {}
    
    // Destructor
    
    ~CscMatrix() {}
    
    CooMatrix tocoo() const;
    
    cArray dotT(const cArrayView b) const;
    
    // Getters
    
    size_t size() const { return i_.size(); }
    size_t rows() const { return m_; }
    size_t cols() const { return n_; }
    lArray const & p() const { return p_; }
    lArray const & i() const { return i_; }
    cArray const & x() const { return x_; }
    
    CscMatrix & operator *= (double r);
    
    CscMatrix & operator &= (const CscMatrix&  B);
    
    CscMatrix& operator ^= (const CscMatrix&  B);
    
    bool hdfsave(const char* name) const;
    
    bool hdfload(const char* name);
    
private:
    
    // dimensions
    long m_;
    long n_;
    
    // representation
    lArray p_;
    lArray i_;
    cArray x_;
    
};

class CsrMatrix {
    
public:
    
    // Constructors
    
    CsrMatrix()
        : m_(0), n_(0), name_() {}
    CsrMatrix(size_t m, size_t n)
        : m_(m), n_(n), name_() {}
    CsrMatrix(CsrMatrix const & A)
        : m_(A.m_), n_(A.n_), p_(A.p_), i_(A.i_), x_(A.x_), name_() {}
    CsrMatrix(size_t m, size_t n, lArrayView const & p, lArrayView const & i, cArrayView const & x)
        : m_(m), n_(n), p_(p), i_(i), x_(x), name_() {}
    
    // Destructor
    
    ~CsrMatrix() {}
    
    void drop ()
    {
        m_ = n_ = 0;
        p_.drop();
        i_.drop();
        x_.drop();
    }
    
    cArray dot(cArrayView const & b) const;
    
    // Getters
    
    size_t size() const { return i_.size(); }
    size_t rows() const { return m_; }
    size_t cols() const { return n_; }
    lArray const & p() const { return p_; }
    lArray const & i() const { return i_; }
    cArray const & x() const { return x_; }
    
    // return absolute value of the (in absolute value) largest element
    double norm() const;
    
    void write(const char* filename) const;
    
#ifndef NO_PNG

    class PngGenerator : public png::generator<png::gray_pixel_1,PngGenerator>
    {
        typedef png::generator<png::gray_pixel_1,PngGenerator> base_t;
        typedef png::packed_pixel_row<png::gray_pixel_1> row;
        typedef png::row_traits<row> row_traits;
        
        public:
            
            PngGenerator (CsrMatrix const * mat, double threshold);
            ~PngGenerator ();
            
            png::byte* get_next_row (size_t pos);
            
        private:
            
            CsrMatrix const * M_;
            row buff_;
            double threshold_;
    };

    void plot (const char* filename, double threshold = 0.) const;
#endif

#ifndef NO_UMFPACK

    class LUft
    {
    public:
        
        LUft ()
            : numeric_(nullptr), matrix_(nullptr), filename_(), info_(UMFPACK_INFO) {}
        
        LUft (LUft const & lu)
            : numeric_(lu.numeric_), matrix_(lu.matrix_), filename_(lu.filename_), info_(lu.info_) {}
        
        LUft (LUft && lu)
            : numeric_(lu.numeric_), matrix_(lu.matrix_), filename_(lu.filename_), info_(lu.info_) { lu.numeric_ = nullptr; }
        
        LUft (const CsrMatrix* matrix, void* numeric)
            : numeric_(numeric), matrix_(matrix), filename_(), info_(UMFPACK_INFO) {}
        
        ~LUft ()
        {
            drop ();
        }
        
        void transfer (LUft && lu)
        {
            // clean this object
            drop();
            
            // transfer data
            numeric_  = lu.numeric_;
            matrix_   = lu.matrix_;
            filename_ = lu.filename_;
            info_     = lu.info_;
            
            // clean source
            lu.numeric_  = nullptr;
        }
        
        size_t size () const
        {
            long lnz, unz, m, n, nz_udiag;
            long status = umfpack_zl_get_lunz (
                &lnz, &unz, &m, &n, &nz_udiag, numeric_
            );
            return status == 0 ? (lnz + unz) * 16 : 0; // Byte count
        }
        
        cArray solve (const cArrayView b, unsigned eqs = 1) const
        {
            // reserve space for the solution
            cArray x(b.size());
            
            // solve
            solve(b, x, eqs);
            
            // return the result
            return x;
        }
        void solve (const cArrayView b, cArrayView x, int eqs = 1) const
        {
            // check sizes
            assert (eqs * matrix_->n_ == (int)x.size());
            assert (eqs * matrix_->n_ == (int)b.size());
            
            // solve for all RHSs
            for (int eq = 0; eq < eqs; eq++)
            {
                // solve for current RHS
                long status = umfpack_zl_solve (
                    UMFPACK_Aat,
                    matrix_->p_.data(), matrix_->i_.data(),
                    reinterpret_cast<const double*>(matrix_->x_.data()), nullptr,
                    reinterpret_cast<double*>(&x[0] + eq * matrix_->n_), nullptr,
                    reinterpret_cast<const double*>(&b[0] + eq * matrix_->n_), nullptr,
                    numeric_, nullptr, &info_[0]
                );
                
                // check output
                if (status != UMFPACK_OK)
                {
                    std::cerr << "\n[CsrMatrix::LUft::solve] Exit status " << status << "\n";
                    umfpack_zl_report_status(0, status);
                }
            }
        }
        
        rArray const & info() const { return info_; }
        
        void link (std::string name) { filename_ = name; }
        
        void save (std::string name) const
        {
            long err = umfpack_zl_save_numeric (numeric_, const_cast<char*>(filename_.c_str()));
            
            if (err == UMFPACK_ERROR_invalid_Numeric_object)
                throw exception ("[LUft::save] Invalid numeric object.");
            
            if (err == UMFPACK_ERROR_file_IO)
                throw exception ("[LUft::save] I/O error.");
        }
        void save () const
        {
            save (filename_);
        }
        
        void load (std::string name)
        {
            long err = umfpack_zl_load_numeric (&numeric_, const_cast<char*>(filename_.c_str()));
            
            if (err == UMFPACK_ERROR_out_of_memory)
                throw exception ("[LUft::load] Out of memory.");
            
            if (err == UMFPACK_ERROR_file_IO)
                throw exception ("[LUft::save] I/O error.");
        }
        void load ()
        {
            load (filename_);
        }
        
        void drop ()
        {
            if (numeric_ != nullptr)
            {
                umfpack_zl_free_numeric (&numeric_);
                numeric_ = nullptr;
            }
        }
        
    private:
        
        void* numeric_;
        
        const CsrMatrix* matrix_;
        
        std::string filename_;
        
        mutable rArray info_;
        
        // Disable bitwise copy
        LUft const & operator= (LUft const &);
    };

    LUft factorize (double droptol = 0) const;
#endif
    
    cArray solve (const cArrayView b, size_t eqs = 1) const;
    
    void link (std::string name);
    
    bool hdfsave () const { return hdfsave(name_); }
    bool hdfsave (std::string name) const;
    
    bool hdfload () { return hdfload(name_); }
    bool hdfload (std::string name);
    
    Complex operator() (unsigned i, unsigned j) const;
    
    CsrMatrix & operator *= (Complex r);
    
    CsrMatrix & operator &= (CsrMatrix const & B);
    
    CsrMatrix& operator ^= (CsrMatrix const & B);
    
    CsrMatrix sparse_like (CsrMatrix const & B) const;
    
    cArray diag() const;
    
    CooMatrix tocoo() const;
    
    RowMatrix<Complex> torow() const;
    
    cArray upperSolve (cArrayView const & b) const;
    
    cArray lowerSolve (cArrayView const & b) const;
    
    template <class Functor> CsrMatrix nzTransform (Functor f) const
    {
        // create output matrix
        CsrMatrix A = *this;
        
        // loop over rows
        for (size_t row = 0; row < rows(); row++)
        {
            // get relevant columns
            size_t idx1 = p_[row];
            size_t idx2 = p_[row + 1];
            
            // loop over columns
            for (size_t idx = idx1; idx < idx2; idx++)
            {
                // which column is this?
                size_t col = i_[idx];
                
                // transform the output matrix
                A.x_[idx] = f(row, col, x_[idx]);
            }
        }
        
        return A;
    }
    
private:
    
    // dimensions
    long m_;
    long n_;
    
    // representation
    lArray p_;
    lArray i_;
    cArray x_;
    
    // linked HDF file
    std::string name_;
};


class CooMatrix {
    
public:
    
    // Empty constructors
    
    CooMatrix()
        : m_(0), n_(0), sorted_(true) {}
    CooMatrix(size_t m, size_t n)
        : m_(m), n_(n), sorted_(true) {}
    CooMatrix(CooMatrix const & A)
        : m_(A.m_), n_(A.n_), i_(A.i_), j_(A.j_), x_(A.x_), sorted_(false) {}
    CooMatrix(size_t m, size_t n, NumberArray<long> const & i, NumberArray<long> const & j, NumberArray<Complex> const & x)
        : m_(m), n_(n), i_(i), j_(j), x_(x), sorted_(false) {}
    
    template <class T> CooMatrix (size_t m, size_t n, T a) : m_(m), n_(n), sorted_(false)
    {
        // initialize from column-major formatted input
        size_t i = 0;
        for (size_t col = 0; col < n; col++)
        {
            for (size_t row = 0; row < m; row++)
            {
                // get element from array
                Complex val = *(a + i);
                
                // if nonzero, store
                if (val != 0.)
                {
                    i_.push_back(row);
                    j_.push_back(col);
                    x_.push_back(val);
                }
                
                // move to next element
                i++;
            }
        }
    }
    
    // Destructor
    ~CooMatrix() {}
    
    operator Complex () const
    {
        if (m_ == 1 and n_ == 1)
        {
            // get number of nonzero matrix elements
            size_t elems = this->shake().x_.size();
            if (elems > 1)
                throw "[CooMatrix::operator Complex] more elements than nominal volume!\n";
            else if (elems == 1)
                return this->shake().x_[0];
            else
                return 0.;
        }
        else
            throw "[CooMatrix::operator Complex] matrix is not 1×1!\n";
    }
    
    // Getters
    
    size_t rows() const { return m_; }
    size_t cols() const { return n_; }
    size_t size() const { return i_.size(); }
    lArray const & i() const { return i_; }
    lArray const & j() const { return j_; }
    cArray const & v() const { return x_; }
    
    Complex operator() (size_t ix, size_t iy) const
    {
        for (size_t n = 0; n < i_.size(); n++)
             if ((size_t)i_[n] == ix and (size_t)j_[n] == iy)
                return x_[n];
        return 0.;
    }
    
    template <class Functor> CooMatrix& symm_populate_band (size_t d, Functor f)
    {
        Complex val;
        
        for (size_t row = 0; row < m_; row++)
        {
            for (size_t col = row; col < n_ and col - row <= d; col++)
            {
                val = f(row,col);
                
                if (val != 0.)
                {
                    i_.push_back(row);
                    j_.push_back(col);
                    x_.push_back(val);
                    
                    if (row != col)
                    {
                        i_.push_back(col);
                        j_.push_back(row);
                        x_.push_back(val);
                    }
                }
            }
        }
        
        sorted_ = false;
        
        return *this;
    }
    
    template <class Functor> CooMatrix& populate (Functor f)
    {
        Complex val;
        
        for (size_t row = 0; row < m_; row++)
        {
            for (size_t col = 0; col < n_; col++)
            {
                val = f(row,col);
                
                if (val != 0.)
                {
                    i_.push_back(row);
                    j_.push_back(col);
                    x_.push_back(val);
                }
            }
        }
        
        sorted_ = false;
        
        return *this;
    }
    
    // Assignment
    CooMatrix & operator = (CooMatrix const & A)
    {
        m_ = A.m_;
        n_ = A.n_;
        i_ = A.i_;
        j_ = A.j_;
        x_ = A.x_;
        
        sorted_ = A.sorted_;
        
        return *this;
    }
    
    void add (long i, long j, Complex v)
    {
        i_.push_back(i);
        j_.push_back(j);
        x_.push_back(v);
        
        sorted_ = false;
    }
    
    CooMatrix transpose() const
    {
        CooMatrix tr;
        
        tr.m_ = n_;
        tr.n_ = m_;
        tr.i_ = j_;
        tr.j_ = i_;
        tr.x_ = x_;
        
        tr.sorted_ = false;
        
        return tr;
    }
    
    CooMatrix& operator += (CooMatrix const & A)
    {
        assert(m_ == A.m_);
        assert(n_ == A.n_);
        
        i_.append(A.i_.begin(), A.i_.end());
        j_.append(A.j_.begin(), A.j_.end());
        x_.append(A.x_.begin(), A.x_.end());
        
        sorted_ = false;
        
        return *this;
    }
    
    CooMatrix& operator -= (CooMatrix const & A)
    {
        assert(m_ == A.m_);
        assert(n_ == A.n_);
        
        size_t prev_size = x_.size();
        
        i_.append(A.i_.begin(), A.i_.end());
        j_.append(A.j_.begin(), A.j_.end());
        x_.append(A.x_.begin(), A.x_.end());
        
        // negate the newly added elements
        for (size_t i = prev_size; i < x_.size(); i++)
            x_[i] = -x_[i];
        
        sorted_ = false;
        
        return *this;
    }
    
    CooMatrix& operator *= (Complex c)
    {
        long nz = i_.size();
        for (long i = 0; i < nz; i++)
            x_[i] *= c;
        
        return *this;
    }
    
    CooMatrix dot (const cArrayView B) const;
    
    Complex ddot (CooMatrix const & B) const;
    
    CooMatrix& operator *= (const cArrayView B);
    
    CooMatrix& operator /= (Complex c)
    {
        return *this *= 1./c;
    }
    
    void resize (size_t m, size_t n)
    {
        m_ = m;
        n_ = n;
    }
    
    CooMatrix reshape (size_t m, size_t n) const;
    
    cArray todense () const;
    
    void sort();
    bool sorted() const { return sorted_; }
    
    CscMatrix tocsc() const;
    
    CsrMatrix tocsr() const;
    
    template <typename DenseMatrixType> DenseMatrixType todense() const
    {
        DenseMatrixType M (rows(), cols());
        for (unsigned idx = 0; idx < x_.size(); idx++)
            M (i_[idx], j_[idx]) = x_[idx];
        return M;
    }
    
    RowMatrix<Complex> torow() const
    {
        return todense<RowMatrix<Complex>>();
    }
    
    ColMatrix<Complex> tocol() const
    {
        return todense<ColMatrix<Complex>>();
    }
    
    SymDiaMatrix todia(MatrixTriangle triangle = lower) const;
    
    cArray solve (const cArrayView b, size_t eqs = 1) const
    {
        // COO format is not optimal for solving -> covert to CSC
        return tocsr().solve(b, eqs);
    }
    
    void write(const char* filename) const;
    
    CooMatrix shake() const;
    
    bool hdfsave(const char* name) const;
    
    bool hdfload(const char* name);

private:

    // dimensions
    size_t m_, n_;

    // ijv-representation
    lArray i_, j_;
    cArray x_;
    
    bool sorted_;
};

class SymDiaMatrix {

public:
    //
    // Constructors
    //

    SymDiaMatrix();
    
    SymDiaMatrix(int n);
    
    SymDiaMatrix(int n, const iArrayView id, const cArrayView v);
    
    SymDiaMatrix(SymDiaMatrix const & A);

    SymDiaMatrix(SymDiaMatrix && A);

    template <class Functor> SymDiaMatrix & populate(unsigned d, Functor f)
    {
        // throw away old data
        elems_.resize(0);
        
        // for all diagonals
        for (size_t id = 0; id <= d; id++)
        {
            // add this diagonal
            idiag_.push_back(id);
            
            // for all elements of the diagonal
            for (int icol = id; icol < n_; icol++)
            {
                // get the row index, too
                int irow = icol - id;
                
                // evaluate the element
                elems_.push_back(f(irow,icol));
            }
        }
        
        return *this;
    }

    //
    // Destructor
    //

    ~SymDiaMatrix() {}
    
    void drop ()
    {
        n_ = 0;
        elems_.drop();
        idiag_.drop();
        dptrs_.clear();
    }
    
    //
    // Getters
    //
    
    iArray const & diag() const { return idiag_; }
    iArray & diag() { return idiag_; }
    
    int diag(int i) const { return idiag_[i]; }
    
    cArrayView main_diagonal() const { return cArrayView(elems_, 0, n_); }
    cArrayView main_diagonal() { return cArrayView(elems_, 0, n_); }
    
    cArray const & data() const { return elems_; }
    cArray & data() { return elems_; }
    
    Complex const * dptr(int i) const { return dptrs_[i]; }
    Complex * dptr(int i) { return dptrs_[i]; }
    
    size_t size() const { return n_; }
    
    int bandwidth() const { return 1 + 2 * idiag_.back(); }
    
    bool is_compatible (SymDiaMatrix const & B) const;
    
    //
    // Arithmetic and other operators
    //
    
    SymDiaMatrix const & operator = (SymDiaMatrix && A);
    SymDiaMatrix const & operator = (SymDiaMatrix const & A);
    
    SymDiaMatrix const & operator += (SymDiaMatrix const & B);
    SymDiaMatrix const & operator -= (SymDiaMatrix const & B);
    
    cArray dot (const cArrayView B, MatrixTriangle triangle = both, bool parallelize = false) const;
    
    cArray lowerSolve (const cArrayView b) const;
    
    cArray upperSolve (const cArrayView b) const;
    
    SymDiaMatrix kron (SymDiaMatrix const & B) const;

    //
    // HDF interface
    //
    
    void link (std::string name);
    
    std::string linkedto () const { return name_; }
    
    bool hdfload () { return hdfload (name_); }
    bool hdfload (std::string name);
    
    bool hdfsave () const { return hdfsave (name_); }
    bool hdfsave (std::string name, bool docompress = false, int consec = 10) const;

    //
    // Conversions to other formats
    //
    
    cArray toPaddedRows () const;
    
    cArray toPaddedCols () const;
    
    CooMatrix tocoo (MatrixTriangle triangle = both) const;
    
    RowMatrix<Complex> torow (MatrixTriangle triangle = both) const;
    
    friend std::ostream & operator << (std::ostream & out, SymDiaMatrix const & A);
    
private:

    // dimension (only square matrices allowed)
    int n_;

    // diagonals: concatenated diagonals starting from the longest
    // to the shortest (i.e. with rising right index)
    cArray elems_;
    
    // upper diagonal indices starting from zero
    iArray idiag_;
    
    // name of linked HDF file
    std::string name_;
    
    // diagonal data pointers
    std::vector<Complex*> dptrs_;
    
    void setup_dptrs_();
};

// --------------------------------------------------------------------------//
// ---- Binary arithmetic operators ---------------------------------------- //
// --------------------------------------------------------------------------//

inline CsrMatrix operator & (const CsrMatrix& A, const CsrMatrix& B)
{
    CsrMatrix C = A;
    return C &= B;
}

inline CsrMatrix operator ^ (const CsrMatrix& A, const CsrMatrix& B)
{
    CsrMatrix C = A;
    return C ^= B;
}

inline CsrMatrix operator * (Complex z, const CsrMatrix&  B)
{
    CsrMatrix C = B;
    return C *= z;
}

inline CsrMatrix operator * (const CsrMatrix& A, double r)
{
    CsrMatrix C = A;
    return C *= r;
}

inline CscMatrix operator & (const CscMatrix& A, const CscMatrix& B)
{
    CscMatrix C = A;
    return C &= B;
}

inline CscMatrix operator ^ (const CscMatrix& A, const CscMatrix& B)
{
    CscMatrix C = A;
    return C ^= B;
}

inline CscMatrix operator * (double r, const CscMatrix& B)
{
    CscMatrix C = B;
    return C *= r;
}

inline CscMatrix operator * (const CscMatrix& A, double r)
{
    CscMatrix C = A;
    return C *= r;
}

inline CooMatrix operator + (const CooMatrix& A, const CooMatrix& B)
{
    CooMatrix C = A;
    return C += B;
}

inline CooMatrix operator - (const CooMatrix& A, const CooMatrix& B)
{
    CooMatrix C = A;
    return C -= B;
}

inline CooMatrix operator * (const Complex& z, const CooMatrix& B)
{
    CooMatrix C = B;
    return C *= z;
}

inline CooMatrix operator * (CooMatrix const & A, const cArrayView B)
{
    CooMatrix C = A;
    return C *= B;
}

inline CooMatrix operator * (const CooMatrix& A, const Complex& z)
{
    CooMatrix C = A;
    return C *= z;
}

SymDiaMatrix operator + (SymDiaMatrix const & A, SymDiaMatrix const & B);
SymDiaMatrix operator - (SymDiaMatrix const & A, SymDiaMatrix const & B);
SymDiaMatrix operator * (SymDiaMatrix const & A, SymDiaMatrix const & B);
SymDiaMatrix operator * (double z, SymDiaMatrix const & A);
SymDiaMatrix operator * (Complex z, SymDiaMatrix const & A);

CooMatrix kron(const CooMatrix& A, const CooMatrix& B);

CooMatrix eye(size_t N);

CooMatrix stairs(size_t N);

#endif