SmithDecomposition.h

Go to the documentation of this file.

/* This file is part of gbLAB.
 *
 * gbLAB is distributed without any warranty under the MIT License.
 */
 
 
#ifndef gbLAB_SmithDecomposition_h_
#define gbLAB_SmithDecomposition_h_
 
#include <utility>      // std::pair, std::make_pair
#include <Eigen/Core>
 
namespace gbLAB
{
    template <int N>
    class SmithDecomposition
    {
        typedef long long int IntValueType;
        typedef Eigen::Matrix<IntValueType,N,N> MatrixNi;
        
        /**********************************************************************/
        static IntValueType gcd(const IntValueType& a,const IntValueType& b)
        {
            return b>0? gcd(b, a % b) : a;
        }
        
        /**********************************************************************/
        template <typename Derived>
        static IntValueType gcd(const Eigen::MatrixBase<Derived>& v)
        {
            IntValueType g=v(0);
            for(int n=1;n<v.size();++n)
            {
                g=gcd(g,v(n));
            }
            return g;
        }
        
        /**********************************************************************/
        static std::pair<int,int> findNonZero(const MatrixNi& A)
        {
            for(int i=0;i<N;++i)
            {
                for(int j=0;j<N;++j)
                {
                    if(A(i,j)!=0)
                    {
                        return std::make_pair(i,j);
                    }
                }
            }
            return std::make_pair(-1,-1);
        }
        
        /**********************************************************************/
        static std::pair<int,int> findNonDivisible(const MatrixNi& A)
        {
            for(int i=1;i<N;++i)
            {
                for(int j=1;j<N;++j)
                {
                    if(A(i,j)%A(0,0)) // D(i,j) is no divisible by D(0,0)
                    {
                        return std::make_pair(i,j);
                    }
                }
            }
            return std::make_pair(-1,-1);
        }
        
        /**********************************************************************/
        void  row_signChange(const int& ID)
        {/* unimodular matrix changing the sign of row/col ID
          */
            MatrixNi T(MatrixNi::Identity());
            T(ID,ID)=-1;
            U=T*U;
            X=X*T;
            D=T*D;
        }
        
        /**********************************************************************/
        void  col_signChange(const int& ID)
        {/* unimodular matrix changing the sign of row/col ID
          */
            MatrixNi T(MatrixNi::Identity());
            T(ID,ID)=-1;
            V=V*T;
            Y=T*Y;
            D=D*T;
        }
        
        /**********************************************************************/
        void row_swap(const int& i, const int& j)
        {/* unimodular matrix swapping rows/cols i and j
          */
            MatrixNi T(MatrixNi::Identity());
            T(i,i)=0;
            T(j,j)=0;
            T(i,j)=1;
            T(j,i)=1;
            U=T*U;
            X=X*T;
            D=T*D;
        }
        
        /**********************************************************************/
        void col_swap(const int& i, const int& j)
        {/* unimodular matrix swapping rows/cols i and j
          */
            MatrixNi T(MatrixNi::Identity());
            T(i,i)=0;
            T(j,j)=0;
            T(i,j)=1;
            T(j,i)=1;
            V=V*T;
            Y=T*Y;
            D=D*T;
        }
        
        /**********************************************************************/
        void row_add_multiply(const int& i, const int& j,const IntValueType& n)
        {/* unimodular matrix adding n times row (col) j to row (col) i
          */
            // note that add applied to colums needs transpose!
            MatrixNi T(MatrixNi::Identity());
            T(i,j)=n;
            U=T*U;
            D=T*D;
            T(i,j)=-n;
            X=X*T;
        }
        
        /**********************************************************************/
        void col_add_multiply(const int& i, const int& j,const IntValueType& n)
        {/* unimodular matrix adding n times row (col) j to row (col) i
          */
            // note that add applied to colums needs transpose!
            MatrixNi T(MatrixNi::Identity());
            T(j,i)=n;
            V=V*T;
            D=D*T;
            T(j,i)=-n;
            Y=T*Y;
        }
        
        /**********************************************************************/
        bool reduceFirstRow()
        {
            // Operate on first row so that D(0,j) is divisible by D(0,0)
            bool didOperate=false;
            int j=1;
            while (j<N)
            {
                const IntValueType r=D(0,j)%D(0,0);
                const IntValueType q=D(0,j)/D(0,0);
                if(r)
                {
                    col_add_multiply(j,0,-q);
                    col_swap(0,j);
                    j=1;
                    didOperate=true;
                }
                else
                {
                    j++;
                }
            }
            return didOperate;
        }
        
        /**********************************************************************/
        bool reduceFirstCol()
        {
            // Operate on first row so that D(0,j) is divisible by D(0,0)
            bool didOperate=false;
            int j=1;
            while (j<N)
            {
                const IntValueType r=D(j,0)%D(0,0);
                const IntValueType q=D(j,0)/D(0,0);
                if(r)
                {
                    row_add_multiply(j,0,-q);
                    row_swap(0,j);
                    j=1;
                    didOperate=true;
                }
                else
                {
                    j++;
                }
            }
            return didOperate;
        }
 
        MatrixNi D;
        MatrixNi U;
        MatrixNi V;
        MatrixNi X;
        MatrixNi Y;
 
        
    public:
                
        /**********************************************************************/
        SmithDecomposition(const MatrixNi& A) : // D=U*A*V, A=X*D*Y
        /* init */ D(A),
        /* init */ U(MatrixNi::Identity()),
        /* init */ V(MatrixNi::Identity()),
        /* init */ X(MatrixNi::Identity()),
        /* init */ Y(MatrixNi::Identity())
        {
            std::pair<int,int> rc=std::make_pair(0,0);
            while (rc.first>=0)
            {
                std::pair<int,int> ij=findNonZero(D);
                if(ij.first>=0)
                {
                    // Swap row(i) and col (j) to amke D(0,0) non-zero.
                    row_swap(0,ij.first);
                    col_swap(0,ij.second);
                    
                    // Operate on first row and col so that all entried are divisible by D(0,0)
                    bool didOperate=true;
                    while (didOperate)
                    {
                        didOperate=reduceFirstRow()+reduceFirstCol();
                    }
                    
                    // If D(0,0) is negative make it positive
                    if(D(0,0)<0)
                    {
                        row_signChange(0);
                    }
                    
                    // Subtract suitable multiples of the first row from the other
                    // rows to replace each entry in the first column, except a(0,0), by zero
                    for (int i=1;i<N;++i)
                    {
                        const IntValueType q=D(i,0)/D(0,0);
                        row_add_multiply(i,0,-q);
                    }
                    
                    // Subtract suitable multiples of the first column from the other
                    // columns to replace each entry in the first row, except a(0,0), by zero
                    for (int j=1;j<N;++j)
                    {
                        const IntValueType q=D(0,j)/D(0,0);
                        col_add_multiply(j,0,-q);
                    }
                    
                    // Now D is in the form [D(0,0) 0s;0s D'] where D' is a (N-1)x(N-1) matrix
                    rc=findNonDivisible(D);
                    if(rc.first>0)
                    {
                        row_add_multiply(0,rc.first,1);
                    }
                    
                    
                    
                    
                }
                else
                {
                    rc=std::make_pair(-1,-1); // end main while loop immediately
                }
                
 
            }
 
            // Recursive iteration
            const SmithDecomposition<N-1> sd(D.template block<N-1,N-1>(1,1));
            D.template block<N-1,N-1>(1,1)=sd.matrixD();
            U.template block<N-1,N>(1,0)=sd.matrixU()*U.template block<N-1,N>(1,0);
            V.template block<1,N-1>(0,1)=V.template block<1,N-1>(0,1)*sd.matrixV();
            V.template block<N-1,N-1>(1,1)=V.template block<N-1,N-1>(1,1)*sd.matrixV();
 
            Y.template block<N-1,N>(1,0)=sd.matrixY()*Y.template block<N-1,N>(1,0);
            X.template block<1,N-1>(0,1)=X.template block<1,N-1>(0,1)*sd.matrixX();
            X.template block<N-1,N-1>(1,1)=X.template block<N-1,N-1>(1,1)*sd.matrixX();
 
            
            // Find a non-zero value of A and make it D(0,0)
            
            const IntValueType UAVD((U*A*V-D).squaredNorm());
            const IntValueType XDYA((X*D*Y-A).squaredNorm());
            
            if (UAVD!=0 || XDYA!=0)
            {
                throw std::runtime_error("Smith decomposition failed\n");
            }
            
        }
        
        /**********************************************************************/
        const MatrixNi& matrixU() const
        {
            return U;
        }
        
        /**********************************************************************/
        const MatrixNi& matrixD() const
        {
            return D;
        }
        
        /**********************************************************************/
        const MatrixNi& matrixV() const
        {
            return V;
        }
        
        /**********************************************************************/
        const MatrixNi& matrixX() const
        {
            return X;
        }
        
        /**********************************************************************/
        const MatrixNi& matrixY() const
        {
            return Y;
        }
    };
    
    /**************************************************************************/
    /**************************************************************************/
    template <>
    class SmithDecomposition<1>
    {
        typedef long long int IntValueType;
        typedef Eigen::Matrix<IntValueType,1,1> MatrixNi;
 
        const MatrixNi D;
        const MatrixNi U;
        const MatrixNi V;
        const MatrixNi X;
        const MatrixNi Y;
        
    public:
        
        /**********************************************************************/
        SmithDecomposition(const MatrixNi& A) :
        /* init */ D(A),
        /* init */ U(MatrixNi::Identity()),
        /* init */ V(MatrixNi::Identity()),
        /* init */ X(MatrixNi::Identity()),
        /* init */ Y(MatrixNi::Identity())
 
        {
        }
        /**********************************************************************/
        const MatrixNi& matrixU() const
        {
            return U;
        }
        
        /**********************************************************************/
        const MatrixNi& matrixD() const
        {
            return D;
        }
        
        /**********************************************************************/
        const MatrixNi& matrixV() const
        {
            return V;
        }
        
        /**********************************************************************/
        const MatrixNi& matrixX() const
        {
            return X;
        }
        
        /**********************************************************************/
        const MatrixNi& matrixY() const
        {
            return Y;
        }
    };
}
#endif