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 oILAB {
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;
        }
    };
    } // namespace oILAB
#endif