C plus plus:Tutorials:Matrices

#ifndef MATRIX2_HPP
#define MATRIX2_HPP
 
#include <cmath>
#include <iostream>
#include <vector>
#include <stdexcept>
#include <algorithm>
#include <functional>
 
template <class T>
class Matrix2
{
public:
    explicit Matrix2(int r = 0, int c = 0, T init_with = T())
     : rows_( r < 1 ? 4 : r), columns_( c < 1 ? 4 : c ),
       data_( rows_ * columns_, init_with ) {}
 
    template <typename U>
    explicit Matrix2(Matrix2<U> const &other)
     : rows_( other.rows() ), columns_( other.columns() ),
       data_( &*other, &*other+(rows_*columns_) ) {}
 
    T &operator()(int row, int column) { return data_[column+(columns_*row)]; }
    T const &operator()(int row, int column) const { return data_[column+(columns_*row)]; }
 
    int rows() const { return rows_; }
    int columns() const { return columns_; }
    operator T*() { return data_.empty() ? 0 : &data_[0]; }
    operator T const*() const { return data_.empty() ? 0 : &data_[0]; }
 
    Matrix2 &operator+=(const Matrix2& other) {
        if ( rows() != other.rows() || columns() != other.columns() ) {
            throw std::invalid_argument("To add matrices they must have the same size");
        }
 
        std::transform( data_.begin(), data_.end(),
                        other.data_.begin(),
                        data_.begin(),
                        std::plus<T>() );
        return *this;
    }
    Matrix2 &operator-=(const Matrix2& other) {
        if ( rows() != other.rows() || columns() != other.columns() ) {
            throw std::invalid_argument("To subtract matrices they must have the same size");
        }
 
        std::transform( data_.begin(), data_.end(),
                        other.data_.begin(),
                        data_.begin(),
                        std::minus<T>() );
        return *this;
    }
    template <class U>
    friend const Matrix2<U> operator*(const Matrix2<U>& one, const Matrix2<U>& two);
    Matrix2 &operator*=(const Matrix2& other) {
        return *this = *this * other;
    }
    Matrix2 &operator*=(T const &scale)
    {
        std::transform( data_.begin(), data_.end(),
                        data_.begin(),
                        std::bind2nd( std::multiplies<T>(), scale ) );
        return *this;
    }
    Matrix2 &operator/=(T const &scale)
    {
        std::transform( data_.begin(), data_.end(),
                        data_.begin(),
                        std::bind2nd( std::divides<T>(), scale ) );
        return *this;
    }
 
    void load_identity() {
        if ( rows() != columns() ) {
            throw std::invalid_argument("Identity matrices only exist for square matrices");
        }
 
        for (int i=0; i<rows(); i++) {
            for (int j=0; j<columns(); j++) {
                operator()( i, j ) = (i==j);
            }
        }
    }
 
    template <typename U>
    Matrix2 &operator=(Matrix2<U> const &other) {
        // note that the assign needs to be done first for exception safety.
        data_.assign( other.data_.begin(), other.data_.end() );
        rows_ = other.rows_;
        columns_ = other.columns_;
        return *this;
    }
 
private:
    // note: changing the order here will break both constructors
    int rows_;
    int columns_;
    std::vector<T> data_;
};
 
template <class T>
const Matrix2<T> operator+(Matrix2<T> one, Matrix2<T> const &two) {
    return one += two;
}
template <class T>
const Matrix2<T> operator-(Matrix2<T> one, Matrix2<T> const &two) {
    return one -= two;
}
 
template <class T>
const Matrix2<T> operator*(const Matrix2<T>& one, const Matrix2<T>& two) {
    if ( one.columns() != two.rows() ) {
        throw std::invalid_argument("Matrices can only be multiplied if the first has"
                                    " exactly as many columns as the second has rows");
    }
 
    Matrix2<T> result( one.rows(), two.columns() );
 
    for (int i=0; i<result.rows(); i++)
        for (int j=0; j<result.columns(); j++)
            for (int a=0; a<one.columns(); a++)
                result(i,j) += one(i,a) * two(a,j);
 
    return result;
}
 
template <class T>
const Matrix2<T> operator*(Matrix2<T> one, T const &factor) {
    return one *= factor;
}
template <class T>
const Matrix2<T> operator*(T const &factor, Matrix2<T> one) {
    return one *= factor;
}
template <class T>
const Matrix2<T> operator/(Matrix2<T> one, T const &factor) {
    return one /= factor;
}
 
template <typename T, typename U>
bool operator==(Matrix2<T> const &one, Matrix2<U> const &two) {
    return (one.rows() == two.rows()
        && one.columns() == two.columns()
        && std::equal( &*one, &*one + one.rows()*one.columns(),
                        &*two ));
}
template <typename T, typename U>
bool operator!=(Matrix2<T> const &one, Matrix2<U> const &two) {
    return !( one == two );
}
 
template<class T>
std::ostream& operator<<(std::ostream& os, const Matrix2<T>& mat) {
    for (int i=0; i< mat.rows(); i++) {
        os << "\n[";
        for (int j=0; j<mat.columns(); j++) {
            os << mat(i,j);
            if ( j+1 < mat.columns() )
                os << "\t";
        }
        os << "]" << std::endl;
    }
 
    return os;
}
 
typedef Matrix2<float> Matrixf;
typedef Matrix2<int> Matrixi;
typedef Matrix2<size_t> Matrixu;
 
#endif // MATRIX2_HPP