MMatrix.h

Go to the documentation of this file.

/*
 *  Copyright (c) 2015, Chris Rogers
 *  All rights reserved.
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are met:
 *  1. Redistributions of source code must retain the above copyright notice,
 *     this list of conditions and the following disclaimer.
 *  2. Redistributions in binary form must reproduce the above copyright notice,
 *     this list of conditions and the following disclaimer in the documentation
 *     and/or other materials provided with the distribution.
 *  3. Neither the name of STFC nor the names of its contributors may be used to
 *     endorse or promote products derived from this software without specific
 *     prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 */
 
#ifndef _SRC_COMMON_CPP_MATHS_MMATRIX_HH_
#define _SRC_COMMON_CPP_MATHS_MMATRIX_HH_
 
#include <iostream>
#include <vector>
 
#include "Utilities/GSLBLAS.h"
#include "Utilities/GSLMatrix.h"
 
#include "Utilities/GeneralOpalException.h"
 
#include "Fields/Interpolation/MVector.h"
 
namespace interpolation {
 
 
    template <class Tmplt>
    class MMatrix {
    public:
        MMatrix();
 
        MMatrix(const MMatrix<Tmplt>& mv);
 
        MMatrix(size_t nrows, size_t ncols, Tmplt* data_beginning);
 
        MMatrix(size_t nrows, size_t ncols, Tmplt value);
 
        MMatrix(size_t nrows, size_t ncols);
 
        static MMatrix<Tmplt> Diagonal(size_t i, Tmplt diag_value, Tmplt off_diag_value);
 
        ~MMatrix();
 
        size_t num_row() const;
 
        size_t num_col() const;
 
        Tmplt trace() const;
 
        Tmplt determinant() const;
 
        MMatrix<Tmplt> inverse() const;
 
        void invert();
 
        MMatrix<Tmplt> T() const;
 
        MVector<m_complex> eigenvalues() const;  // return vector of eigenvalues
        std::pair<MVector<m_complex>, MMatrix<m_complex> > eigenvectors()
                const;  // return pair of eigenvalues, eigenvectors (access using my_pair.first or
                        // my_pair.second)
 
        // return a submatrix, with data from min_row to max_row and min_row to max_row inclusive;
        // indexing starts at 1
        MMatrix<Tmplt> sub(size_t min_row, size_t max_row, size_t min_col, size_t max_col) const;
        // create a column vector from the column^th column
        MVector<Tmplt> get_mvector(size_t column) const;
        // return a reference to the datum; indexing starts at 1, goes to num_row (inclusive)
        const Tmplt& operator()(size_t row, size_t column) const;
        Tmplt& operator()(size_t row, size_t column);
 
        MMatrix<Tmplt>& operator=(const MMatrix<Tmplt>& mm);
 
        // TODO - implement iterator
        // class iterator
        //{
        // }
 
        friend MMatrix<m_complex>& operator*=(MMatrix<m_complex>& m, m_complex c);
        friend MMatrix<double>& operator*=(MMatrix<double>& m, double d);
        friend MMatrix<m_complex>& operator*=(MMatrix<m_complex>& m1, MMatrix<m_complex> m2);
        friend MMatrix<double>& operator*=(MMatrix<double>& m1, MMatrix<double> m2);
        friend MVector<m_complex> operator*(MMatrix<m_complex> m, MVector<m_complex> v);
        friend MVector<double> operator*(MMatrix<double> m, MVector<double> v);
        friend MMatrix<m_complex>& operator+=(MMatrix<m_complex>& m1, const MMatrix<m_complex>& m2);
        friend MMatrix<double>& operator+=(MMatrix<double>& m1, const MMatrix<double>& m2);
        template <class Tmplt2>
        friend MMatrix<Tmplt2> operator+(MMatrix<Tmplt2> m1, const MMatrix<Tmplt2> m2);
 
        friend class MMatrix<double>;  // To do the eigenvector problem, MMatrix<double> needs to
                                       // see MMatrix<complex>'s _matrix
 
    private:
        // build the matrix with size i,j, data not initialised
        void build_matrix(size_t i, size_t j);
        // build the matrix with size i,j, data initialised to data in temp
        void build_matrix(size_t i, size_t j, Tmplt* temp);
        // delete the matrix and set it to null
        void delete_matrix();
        // return the matrix overloaded to the correct type or throw
        // if _matrix == nullptr
        static gsl_matrix* get_matrix(const MMatrix<double>& m);
        static gsl_matrix_complex* get_matrix(const MMatrix<m_complex>& m);
        // gsl object
        void* _matrix;
    };
 
    // Operators do what you would expect - i.e. normal mathematical functions
    MMatrix<m_complex>& operator*=(MMatrix<m_complex>& m, m_complex c);
    MMatrix<double>& operator*=(MMatrix<double>& m, double d);
    MMatrix<m_complex>& operator*=(
            MMatrix<m_complex>& m1, MMatrix<m_complex> m2);  // M1 *= M2 returns M1 = M1*M2
    MMatrix<double>& operator*=(
            MMatrix<double>& m1, MMatrix<double> m2);  // M1 *= M2 returns M1 = M1*M2
    MVector<m_complex> operator*(MMatrix<m_complex> m, MVector<m_complex> v);
    MVector<double> operator*(MMatrix<double> m, MVector<double> v);
    MMatrix<m_complex>& operator+=(MMatrix<m_complex>& m1, const MMatrix<m_complex>& m2);
    MMatrix<double>& operator+=(MMatrix<double>& m1, const MMatrix<double>& m2);
 
    template <class Tmplt>
    MMatrix<Tmplt>& operator-=(MMatrix<double>& m1, MMatrix<double> m2);
 
    template <class Tmplt>
    MMatrix<Tmplt> operator*(MMatrix<Tmplt>, MMatrix<Tmplt>);
    template <class Tmplt>
    MMatrix<Tmplt> operator*(MMatrix<Tmplt>, MVector<Tmplt>);
    template <class Tmplt>
    MMatrix<Tmplt> operator*(MMatrix<Tmplt>, Tmplt);
    template <class Tmplt>
    MMatrix<Tmplt> operator/(MMatrix<Tmplt>, Tmplt);
 
    template <class Tmplt>
    MMatrix<Tmplt> operator+(MMatrix<Tmplt>, MMatrix<Tmplt>);
    template <class Tmplt>
    MMatrix<Tmplt> operator-(MMatrix<Tmplt>, MMatrix<Tmplt>);
    template <class Tmplt>
    MMatrix<Tmplt> operator-(MMatrix<Tmplt>);  // unary minus returns matrix*(-1)
 
    template <class Tmplt>
    bool operator==(const MMatrix<Tmplt>& c1, const MMatrix<Tmplt>& c2);
    template <class Tmplt>
    bool operator!=(const MMatrix<Tmplt>& c1, const MMatrix<Tmplt>& c2);
 
    // format is
    //  num_row  num_col
    //  m11 m12 m13 ...
    //  m21 m22 m23 ...
    //  ...
    template <class Tmplt>
    std::ostream& operator<<(std::ostream& out, MMatrix<Tmplt> mat);
    template <class Tmplt>
    std::istream& operator>>(std::istream& in, MMatrix<Tmplt>& mat);
 
    // return matrix of doubles filled with either real or imaginary part of m
    MMatrix<double> re(MMatrix<m_complex> m);
    MMatrix<double> im(MMatrix<m_complex> m);
    // return matrix of m_complex filled with real and imaginary parts
    MMatrix<m_complex> complex(MMatrix<double> real);
    MMatrix<m_complex> complex(MMatrix<double> real, MMatrix<double> imaginary);
 
 
 
    template <>
    inline size_t MMatrix<double>::num_row() const {
        if (_matrix) return ((gsl_matrix*)_matrix)->size1;
        return 0;
    }
    template <>
    inline size_t MMatrix<m_complex>::num_row() const {
        if (_matrix) return ((gsl_matrix_complex*)_matrix)->size1;
        return 0;
    }
    template <>
    inline size_t MMatrix<double>::num_col() const {
        if (_matrix) return ((gsl_matrix*)_matrix)->size2;
        return 0;
    }
 
    template <>
    inline size_t MMatrix<m_complex>::num_col() const {
        if (_matrix) return ((gsl_matrix_complex*)_matrix)->size2;
        return 0;
    }
 
    MMatrix<m_complex> inline& operator*=(MMatrix<m_complex>& m, m_complex c) {
        gsl_matrix_complex_scale((gsl_matrix_complex*)m._matrix, c);
        return m;
    }
 
    MMatrix<double> inline& operator*=(MMatrix<double>& m, double d) {
        gsl_matrix_scale((gsl_matrix*)m._matrix, d);
        return m;
    }
 
    MVector<m_complex> inline operator*(MMatrix<m_complex> m, MVector<m_complex> v) {
        MVector<m_complex> v0(m.num_row());
        gsl_blas_zgemv(
                CblasNoTrans, m_complex_build(1.), (gsl_matrix_complex*)m._matrix,
                (gsl_vector_complex*)v._vector, m_complex_build(0.),
                (gsl_vector_complex*)v0._vector);
        return v0;
    }
 
    MVector<double> inline operator*(MMatrix<double> m, MVector<double> v) {
        MVector<double> v0(m.num_row());
        gsl_blas_dgemv(
                CblasNoTrans, 1., (gsl_matrix*)m._matrix, (gsl_vector*)v._vector, 0.,
                (gsl_vector*)v0._vector);
        return v0;
    }
 
    MMatrix<m_complex> inline& operator+=(MMatrix<m_complex>& m1, const MMatrix<m_complex>& m2) {
        gsl_matrix_complex_add((gsl_matrix_complex*)m1._matrix, (gsl_matrix_complex*)m2._matrix);
        return m1;
    }
 
    MMatrix<double> inline& operator+=(MMatrix<double>& m1, const MMatrix<double>& m2) {
        gsl_matrix_add((gsl_matrix*)m1._matrix, (gsl_matrix*)m2._matrix);
        return m1;
    }
 
    template <class Tmplt>
    MMatrix<Tmplt> inline operator-(MMatrix<Tmplt> m1) {
        for (size_t i = 1; i <= m1.num_row(); i++)
            for (size_t j = 1; j <= m1.num_col(); j++)
                m1(i, j) = -1. * m1(i, j);
        return m1;
    }
 
    template <class Tmplt>
    MMatrix<Tmplt> inline operator*(MMatrix<Tmplt> m1, MMatrix<Tmplt> m2) {
        return m1 *= m2;
    }
 
    template <class Tmplt>
    MMatrix<Tmplt> inline operator*(MMatrix<Tmplt> m, MVector<Tmplt> v) {
        return m *= v;
    }
 
    template <class Tmplt>
    MMatrix<Tmplt> inline operator*(MMatrix<Tmplt> m, Tmplt t) {
        return m *= t;
    }
 
    template <class Tmplt>
    MMatrix<Tmplt> inline& operator/=(MMatrix<Tmplt>& m, Tmplt t) {
        return m *= 1. / t;
    }
    template <class Tmplt>
    MMatrix<Tmplt> inline operator/(MMatrix<Tmplt> m, Tmplt t) {
        return m /= t;
    }
 
    template <class Tmplt>
    MMatrix<Tmplt> inline operator+(MMatrix<Tmplt> m1, MMatrix<Tmplt> m2) {
        MMatrix<Tmplt> m3 = m1 += m2;
        return m3;
    }
 
    template <class Tmplt>
    MMatrix<Tmplt> inline& operator-=(MMatrix<Tmplt>& m1, MMatrix<Tmplt> m2) {
        return m1 += -m2;
    }
    template <class Tmplt>
    MMatrix<Tmplt> inline operator-(MMatrix<Tmplt> m1, MMatrix<Tmplt> m2) {
        return m1 -= m2;
    }
 
    template <>
    const double inline& MMatrix<double>::operator()(size_t i, size_t j) const {
        return *gsl_matrix_ptr((gsl_matrix*)_matrix, i - 1, j - 1);
    }
    template <>
    const m_complex inline& MMatrix<m_complex>::operator()(size_t i, size_t j) const {
        return *gsl_matrix_complex_ptr((gsl_matrix_complex*)_matrix, i - 1, j - 1);
    }
 
    template <>
    double inline& MMatrix<double>::operator()(size_t i, size_t j) {
        return *gsl_matrix_ptr((gsl_matrix*)_matrix, i - 1, j - 1);
    }
    template <>
    m_complex inline& MMatrix<m_complex>::operator()(size_t i, size_t j) {
        return *gsl_matrix_complex_ptr((gsl_matrix_complex*)_matrix, i - 1, j - 1);
    }
 
    template <class Tmplt>
    bool operator==(const MMatrix<Tmplt>& lhs, const MMatrix<Tmplt>& rhs) {
        if (lhs.num_row() != rhs.num_row() || lhs.num_col() != rhs.num_col()) return false;
        for (size_t i = 1; i <= lhs.num_row(); i++)
            for (size_t j = 1; j <= lhs.num_col(); j++)
                if (lhs(i, j) != rhs(i, j)) return false;
        return true;
    }
 
    template <class Tmplt>
    bool inline operator!=(const MMatrix<Tmplt>& lhs, const MMatrix<Tmplt>& rhs) {
        return !(lhs == rhs);
    }
 
    // iostream
    template <class Tmplt>
    std::ostream& operator<<(std::ostream& out, MMatrix<Tmplt>);
    // template <class Tmplt> std::istream& operator>>(std::istream& in,  MMatrix<Tmplt>);
 
    template <class Tmplt>
    gsl_matrix inline* MMatrix<Tmplt>::get_matrix(const MMatrix<double>& m) {
        if (m._matrix == nullptr)
            throw(GeneralOpalException(
                    "MMatrix::get_matrix", "Attempt to access uninitialised matrix"));
        return (gsl_matrix*)m._matrix;
    }
 
    template <class Tmplt>
    gsl_matrix_complex inline* MMatrix<Tmplt>::get_matrix(const MMatrix<m_complex>& m) {
        if (m._matrix == nullptr)
            throw(GeneralOpalException(
                    "MMatrix::get_matrix", "Attempt to access uninitialised matrix"));
        return (gsl_matrix_complex*)m._matrix;
    }
 
}  // namespace interpolation
 
#endif