fixed_symmetric_matrix.h

/* Copyright (C) 2003-2005  Martin Förg
 
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * 
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 *
 * Contact:
 *   martin.o.foerg@googlemail.com
 *
 */
 
#ifndef fixed_symmetric_matrix_h
#define fixed_symmetric_matrix_h
 
#include "types.h"
#include <vector>
 
namespace fmatvec {
 
  template <int M, class AT> class Matrix<Symmetric,Fixed<M>,Fixed<M>,AT> {
 
    protected:
 
 
      AT ele[M][M];
 
      template <class Type, class Row, class Col> inline Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>& copy(const Matrix<Type,Row,Col,AT> &A);
      template <class Row> inline Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>& copy(const Matrix<Symmetric,Row,Row,AT> &A);
 
 
    public:
      static constexpr bool isVector {false};
      using value_type = AT;
      using shape_type = Symmetric;
 
      explicit Matrix(Noinit ini) { }
      explicit Matrix(Init ini=INIT, const AT &a=AT()) { init(a); }
      explicit Matrix(Eye ini, const AT &a=1) { init(ini,a); }
      explicit Matrix(int m, int n, Noinit ini) { FMATVEC_ASSERT(m==M && n==M, AT); }
      explicit Matrix(int m, int n, Init ini, const AT &a=AT()) { FMATVEC_ASSERT(m==M && n==M, AT); init(a); }
      explicit Matrix(int m, int n, Eye ini, const AT &a=1) { FMATVEC_ASSERT(m==M && n==M, AT); init(ini,a); }
 
      Matrix(const Matrix<Symmetric,Fixed<M>,Fixed<M>,AT> &A) = default;
 
      template<class Row>
      Matrix(const Matrix<Symmetric,Row,Row,AT> &A) {
    FMATVEC_ASSERT(A.size() == M, AT); 
    copy(A);
      }
 
      template<class Type, class Row, class Col>
      explicit Matrix(const Matrix<Type,Row,Col,AT> &A) {
    FMATVEC_ASSERT(A.rows() == M, AT); 
    FMATVEC_ASSERT(A.cols() == M, AT);
    copy(A);
      }
 
      inline Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>& operator=(const Matrix<Symmetric,Fixed<M>,Fixed<M>,AT> &A) = default;
 
      template<class Type, class Row, class Col>
      inline Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>& operator=(const Matrix<Type,Row,Col,AT> &A) {
        FMATVEC_ASSERT(A.rows() == A.cols(), AT);
        FMATVEC_ASSERT(M == A.rows(), AT);
        return copy(A);
      }
 
      template<class Type, class Row, class Col>
      inline Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>& operator<<=(const Matrix<Type,Row,Col,AT> &A) { return operator=(A); }
 
      void resize(int n, int m) {
        if(n!=M || m!=M)
          throw std::runtime_error("A fixed symmetric matrix cannot be resized.");
      }
 
      AT& operator()(int i, int j) {
    FMATVEC_ASSERT(i>=0, AT);
    FMATVEC_ASSERT(j>=0, AT);
    FMATVEC_ASSERT(i<M, AT);
    FMATVEC_ASSERT(j<M, AT);
 
    return e(i,j);
      }
 
      const AT& operator()(int i, int j) const {
    FMATVEC_ASSERT(i>=0, AT);
    FMATVEC_ASSERT(j>=0, AT);
    FMATVEC_ASSERT(i<M, AT);
    FMATVEC_ASSERT(j<M, AT);
 
    return e(i,j);
      }
 
      AT& ei(int i, int j) {
    return ele[i][j];
      }
 
      const AT& ei(int i, int j) const {
    return ele[i][j];
      }
 
      AT& ej(int i, int j) {
    return ele[j][i];
      }
 
      const AT& ej(int i, int j) const {
    return ele[j][i];
      }
 
       AT& e(int i, int j) {
    return j > i ? ej(i,j) : ei(i,j);
      }
 
      const AT& e(int i, int j) const {
    return j > i ? ej(i,j) : ei(i,j);
      }
 
      AT* operator()() {return ele;}
 
      const AT* operator()() const {return ele;}
 
      constexpr int size() const {return M;}
 
      constexpr int rows() const {return M;}
 
      constexpr int cols() const {return M;}
 
      int ldim() const {return M;}
 
      CBLAS_ORDER blasOrder() const {
        return CblasRowMajor;
      }
 
      CBLAS_UPLO blasUplo() const {
        return  CblasLower;
      }
 
      inline Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>& init(const AT &val=AT()); 
      inline Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>& init(Init, const AT &a=AT()) { return init(a); }
      inline Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>& init(Eye, const AT &val=1);
      inline Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>& init(Noinit, const AT &a=AT()) { return *this; }
 
      explicit inline operator std::vector<std::vector<AT>>() const;
 
      explicit inline Matrix(const std::vector<std::vector<AT>> &m);
 
//      /*! \brief Cast to AT.
//       *
//       * \return The AT representation of the matrix
//       * */
//      explicit operator AT() const {
//        FMATVEC_ASSERT(M==1, AT);
//        return ele[0][0];
//      }
//
//      /*! \brief AT Constructor.
//       * Constructs and initializes a matrix with a AT object.
//       * \param x The AT the matrix will be initialized with.
//       * */
//      explicit Matrix(const AT &x) {
//        FMATVEC_ASSERT(M==1, AT);
//        ele[0][0] = x;
//      }
  };
 
  template <int M, class AT>
    inline Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>& Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>::init(const AT &val) {
 
      for(int i=0; i<M; i++) 
        for(int j=i; j<M; j++) 
          ej(i,j) = val;
 
      return *this;
    }
 
  template <int M, class AT>
    inline Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>& Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>::init(Eye, const AT &val) {
 
      for(int i=0; i<size(); i++) {
        ej(i,i) = val;
        for(int j=i+1; j<size(); j++) {
          ej(i,j) = 0; 
        }
      }
      return *this;
    }
 
  template <int M, class AT>
    inline Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>::operator std::vector<std::vector<AT>>() const {
      std::vector<std::vector<AT>> ret(rows(),std::vector<AT>(cols()));
      for(int r=0; r<rows(); r++) {
        for(int c=r; c<cols(); c++) {
          ret[r][c]=ej(r,c);
      if(c>r)
        ret[c][r]=ej(r,c);
    }
      }
      return ret;
    }
 
  template <int M, class AT>
    inline Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>::Matrix(const std::vector<std::vector<AT>> &m) {
      if(static_cast<int>(m.size())!=rows())
        throw std::runtime_error("The number of rows does not match.");
      for(int r=0; r<rows(); r++) {
        if(static_cast<int>(m[r].size())!=cols())
          throw std::runtime_error("The rows of the input have different length.");
        for(int c=r; c<cols(); c++) {
          using std::abs;
          ej(r,c)=m[r][c];
          if(c>r && abs(m[r][c]-m[c][r])>abs(m[r][c]*1e-13+1e-13))
            throw std::runtime_error("The input is not symmetric.");
        }
      }
    }
 
 
  template <int M, class AT> template <class Type, class Row, class Col>
    inline Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>& Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>::copy(const Matrix<Type,Row,Col,AT> &A) {
      for(int i=0; i<M; i++) 
        for(int j=i; j<M; j++)
          ej(i,j) = A.e(i,j);
      return *this;
    }
 
  template <int M, class AT> template <class Row>
    inline Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>& Matrix<Symmetric,Fixed<M>,Fixed<M>,AT>::copy(const Matrix<Symmetric,Row,Row,AT> &A) {
      for(int i=0; i<M; i++) 
        for(int j=i; j<M; j++)
          ej(i,j) = A.ej(i,j);
      return *this;
    }
 
 
}
 
#endif