libigl_Martin/include/igl/kronecker_product.cpp

// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2014 Daniele Panozzo <daniele.panozzo@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla Public License
// v. 2.0. If a copy of the MPL was not distributed with this file, You can
// obtain one at http://mozilla.org/MPL/2.0/.
#include "kronecker_product.h"

//// Bug in unsupported/Eigen/SparseExtra needs iostream first
//#include <iostream>
//#include <unsupported/Eigen/SparseExtra>
//
//template <typename Scalar>
//IGL_INLINE Eigen::SparseMatrix<Scalar> igl::kronecker_product(
//  const Eigen::SparseMatrix<Scalar> & A,
//  const Eigen::SparseMatrix<Scalar> & B)
//{
//  using namespace Eigen;
//  using namespace std;
//
//  // Convert B in triplets format
//  MatrixXd B_triplets(B.nonZeros(),3);
//  int count = 0;
//  for (int k=0; k<B.outerSize(); ++k)
//    for (SparseMatrix<double>::InnerIterator it(B,k); it; ++it)
//      B_triplets.row(count++) << it.row(), it.col(), it.value();
//
//  MatrixXd C_triplets(B_triplets.rows()*A.nonZeros(),3);
//  count = 0;
//  for (int k=0; k<A.outerSize(); ++k)
//    for (SparseMatrix<double>::InnerIterator it(A,k); it; ++it)
//    {
//      int i = it.row();
//      int j = it.col();
//      double v = it.value();
//
//      MatrixXd B_triplets_copy = B_triplets;
//      B_triplets_copy.col(0) = B_triplets_copy.col(0).array() + double(B.rows()*i);
//      B_triplets_copy.col(1) = B_triplets_copy.col(1).array() + double(B.cols()*j);
//      B_triplets_copy.col(2) = B_triplets_copy.col(2).array() * v;
//
//      C_triplets.block(count*B_triplets.rows(),0,
//                       B_triplets.rows(), B_triplets.cols()) = B_triplets_copy;
//
//      count++;
//    }
//
//  typedef Eigen::Triplet<double> T;
//  std::vector<T> triplets;
//  triplets.reserve(C_triplets.rows());
//  
//  for(unsigned i=0; i<C_triplets.rows(); ++i)
//    triplets.push_back(T(C_triplets(i,0),C_triplets(i,1),C_triplets(i,2)));
//  SparseMatrix<Scalar> C(A.rows()*B.rows(),A.cols()*B.cols());
//  C.setFromTriplets(triplets.begin(),triplets.end());
//  
//  return C;
//}
//
//#ifdef IGL_STATIC_LIBRARY
//// Explicit template specialization
//// generated by autoexplicit.sh
//template Eigen::SparseMatrix<double, 0, int> igl::kronecker_product<double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int> const&);
//#endif