SymMatrix.h

// ===========================================================================
// Imagine++ Libraries
// Copyright (C) Imagine
// For detailed information: http://imagine.enpc.fr/software
// ===========================================================================

namespace Imagine {

    template <typename T> class Matrix;

    template <typename T>
    class SymMatrix : public Array<T> {
        typedef Array<T> Base;
        int _n;
        size_t sz(int n) { return size_t(n)*(n+1)/2; } 
    public:
        typedef typename Base::iterator iterator;
        typedef typename Base::const_iterator const_iterator;
        SymMatrix() : Base(), _n(0) {}
        explicit SymMatrix(int N) : Base(sz(N)), _n(N) {}
        SymMatrix(const SymMatrix& A) : Base(A), _n(A._n) {}
        SymMatrix(T* t, int N,bool handleDelete=false) : Base(t,sz(N),handleDelete), _n(N) {  }
        SymMatrix(const Matrix<T>& A) : Base(sz(A.nrow())), _n(A.nrow()) { // Upper part
            assert(A.nrow()==A.ncol());
            for (int j=0; j<A.nrow(); j++)
                equalize(j+1,A.data()+size_t(j)*A.nrow(),1,this->data()+sz(j),1);
        }
        ~SymMatrix() {}
        SymMatrix& operator=(const SymMatrix& A) { Base::operator =(A); _n = A._n; return *this; }
        SymMatrix clone() const
        {
            SymMatrix A(nrow());
            equalize((int)this->size(),this->data(),1,A.data(),1);
            return A;
        }
        void setSize(int n) {
            Base::setSize(sz(n));
            _n = n;
        }
        int nrow() const { return _n; }
        int ncol() const { return _n; }

        SymMatrix& fill(T x) { Base::fill(x); return *this; }
        static SymMatrix Zero(int n){
            SymMatrix Z(n);
            memset(Z.data(),0,Z.size()*sizeof(T));
            return Z;
        }
        static SymMatrix Identity(int n){
            SymMatrix I=SymMatrix::Zero(n);
            for (int i=0;i<n;i++)
                I(i,i)=T(1);
            return I;
        }
        const T& operator()(int i,int j) const {
            assert(i>=0 && i<_n && j>=0 && j<_n);
            if(i<=j)
                return (*this)[i+size_t(j)*(j+1)/2];
            else
                return (*this)[j+size_t(i)*(i+1)/2];
        }
        T& operator()(int i,int j) {
            assert(i>=0 && i<_n && j>=0 && j<_n);
            if(i<=j)
                return (*this)[i+size_t(j)*(j+1)/2];
            else
                return (*this)[j+size_t(i)*(i+1)/2];
        }
        Matrix<T> getSubMat(int i0, int m, int j0, int n) const { 
            assert (i0>=0 && i0+m<=nrow() && j0>=0 && j0+n<=ncol());
            Matrix<T> retMat(m,n);
            for(int j=0;j<n;j++)
                for(int i=0;i<m;i++)
                    retMat(i,j)=(*this)(i0+i,j0+j);
            return retMat;
        }
        Vector<T> getDiagonal() const {
            Vector<T> v(_n);
            for (int i=0;i<_n;i++)
                v[i]=(*this)(i,i);
            return v;
        }
        Vector<T> operator *(const Vector<T> &v) const {
            assert((size_t)_n==v.size());
            Vector<T> y(_n);
            y.fill(T(0));
            int _nn=_n;
            vectorProduct(_nn,T(1),this->data(),v.data(),T(0),y.data());
            return y;
        }
        friend SymMatrix inverse(const SymMatrix& A) {

            SymMatrix invA=A.clone();
            T *pivots=new T[A.nrow()*A.nrow()];
            int info,n=A.nrow();
            // LU Inverse

            symInverse(&n,invA.data(),pivots,&info);
      Matrix<T> inv(pivots,A.nrow(),A.nrow());
      SymMatrix<T> invA2(inv);
            delete[] pivots;
            if(info!=0) {
                std::cerr << "Cannot invert symmetric matrix" << std::endl;
                invA.fill(T(0));
                return invA;
            }
            return invA2;
        }
        friend SymMatrix posDefInverse(const SymMatrix& A) {

            SymMatrix invA=A.clone();
            T *pivots=new T[A.nrow()*A.nrow()];
            int info,n=A.nrow();
            // LU Inverse

            symDPInverse(&n,invA.data(),pivots,&info);
      Matrix<T> inv(pivots,A.nrow(),A.nrow());
      SymMatrix<T> invA2(inv);
            delete[] pivots;
            if(info!=0) {
                std::cerr << "Cannot invert symmetric matrix" << std::endl;
                invA.fill(T(0));
                return invA;
            }
            return invA2;
        }
        friend T det(const SymMatrix& A) 
        {
            int info,n=A.nrow();
      SymMatrix symtemp = A.clone();
            T d =1;
      d = symDeterminant(&n,symtemp.data(),&info);
            if(info!=0) {
                std::cerr << "Cannot compute determinant" << std::endl;
                return 0;
            };
            return d;
        }
        friend void EVD(const SymMatrix& A,Matrix<T>& Q, Vector<T>& Lambda) {
            SymMatrix symtemp = A.clone();
            Lambda = Vector<T>(A.nrow());
            Q = Matrix<T>(A.nrow(),A.nrow());
            int info,n=A.nrow();
            symEigenValues(&n,symtemp.data(),Lambda.data(),Q.data(),&info);
            if(info) {
                std::cerr << "Cannot compute EVD decomposition" << std::endl;
                return;
            }
        }
        Matrix<T> operator *(const SymMatrix& B) const {
            // No blas function
            assert(ncol()==B.nrow());
            Matrix<T> C(nrow(),B.ncol());
            for (int j=0;j<B.ncol();j++) 
                for (int i=0;i<nrow();i++) {
                    C(i,j)=0;
                    for (int k=0;k<ncol();k++)
                        C(i,j)+=(*this)(i,k)*B(k,j);
                }
            return C;
        }
        Matrix<T> operator *(const Matrix<T>& B) const {
            // No blas function
            assert(ncol()==B.nrow());
            Matrix<T> C(nrow(),B.ncol());
            for (int j=0;j<B.ncol();j++) 
                for (int i=0;i<nrow();i++) {
                    C(i,j)=0;
                    for (int k=0;k<ncol();k++)
                        C(i,j)+=(*this)(i,k)*B(k,j);
                }
            return C;
        }
        SymMatrix operator +(const SymMatrix& B) const {
            SymMatrix C=clone();
            return (C+=B);
        }
        SymMatrix& operator +=(const SymMatrix& B) {
            assert(nrow()==B.nrow());
            combine((int)this->size(), T(1), B.data(), this->data());
            return *this;
        }
        SymMatrix operator -(const SymMatrix& B) const {
            SymMatrix C=clone();
            return (C-=B);
        }
        SymMatrix& operator -=(const SymMatrix& B) {
            assert(nrow()==B.nrow());
            combine((int)this->size(), T(-1), B.data(), this->data());
            return *this;
        }
        SymMatrix operator - () const {
            return (*this)*T(-1);
        }
        SymMatrix operator *(T x) const {
            SymMatrix C=clone();
            return (C*=x);
        }
        friend inline SymMatrix operator *(T x, const SymMatrix& B) {return B*x;}
        SymMatrix operator /(T x) const {
            SymMatrix C=clone();
            return (C/=x);
        }
        SymMatrix& operator *=(T x) {
            multiply((int)this->size(), x, this->data());
            return *this;
        }
        SymMatrix& operator /=(T x) {
            multiply((int)this->size(), T(1)/x, this->data());
            return *this;
        }
        friend Vector<T> linSolve(const SymMatrix& A,const Vector<T>& b) 
        {
            assert(b.size()==(size_t)A.nrow());
            SymMatrix invA=A.clone();
            Vector<T> x=b.clone();
            // Bunch Kaufmann factorization
            int *pivots=new int[A.nrow()];
            int info,n=A.nrow();
            // Solve linear system using factorization
            symSolve(&n,invA.data(),pivots,x.data(),&info);
            if (info) {
                x.fill(T(0));
                std::cerr<<"Cannot solve linear system"<<std::endl;
            }
            delete[] pivots;
            return x;
        }
        void write(std::ostream& out) const {
            int n=nrow();
            out.write((const char*)&n,sizeof(int));
            out.write((const char*)this->data(),(std::streamsize)this->size()*sizeof(T));
        }
        void read(std::istream& in) {
            int n;
            in.read((char*)&n,sizeof(int));
            setSize(n);
            in.read((char*)this->data(),(std::streamsize)this->size()*sizeof(T));
        }
        friend inline std::ostream& operator<<(std::ostream& out,const SymMatrix& A) {
            out<<A.nrow()<<" "<<A.ncol()<< std::endl;
            for (int i=0;i<A.nrow();i++) {
                for (int j=0;j<A.ncol();j++) {
                    out<<A(i,j);
                    if (j<A.ncol()-1) out<<" ";
                }
                out<<std::endl;
            }
            return out;
        }
        friend inline std::istream& operator>>(std::istream& in,SymMatrix& A) {
            int m,n;
            in>>m>>n;
            assert(m==n);
            A.setSize(n);
            for (int i=0;i<A.nrow();i++)
                for (int j=0;j<A.ncol();j++)
                    in>>A(i,j);
            return in;
        }
    };

#ifndef DOXYGEN_SHOULD_SKIP_THIS
    // Declaration in Matrix.h
    template <typename T> 
    Matrix<T>::Matrix(const SymMatrix<T>& A) : MultiArray<T,2>(A.nrow(),A.ncol()) {
        for (int j=0; j<ncol(); j++)
            for (int i=0; i<nrow(); i++)
                (*this)(i,j)=A(i,j);
    }
    // Declaration in Matrix.h
    template <typename T> 
    Matrix<T> Matrix<T>::operator *(const SymMatrix<T>& B) const
    {
        assert(ncol()==B.nrow());
        Matrix<T> C(nrow(),B.ncol());
        for (int j=0;j<B.ncol();j++) 
            for (int i=0;i<nrow();i++) {
                C(i,j)=0;
                for (int k=0;k<ncol();k++)
                    C(i,j)+=(*this)(i,k)*B(k,j);
            }
        return C;
    }
#endif

}