blkdiag function

include/igl/blkdiag.cpp

@@ -0,0 +1,77 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2020 Alec Jacobson <[email protected]>
+// 
+// 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 "blkdiag.h"
+
+template <typename Scalar>
+IGL_INLINE void igl::blkdiag(
+  const std::vector<Eigen::SparseMatrix<Scalar>> & L, 
+  Eigen::SparseMatrix<Scalar> & Y)
+{
+  int nr = 0;
+  int nc = 0;
+  int nnz = 0;
+  for(const auto & A : L)
+  {
+    nr += A.rows();
+    nc += A.cols();
+  }
+  Y.resize(nr,nc);
+  Eigen::VectorXi per_col = Eigen::VectorXi::Zero(nc);
+  for(int j = 0;j<L.size();j++)
+  {
+    per_col(j) = L[j].cols();
+  }
+  Y.reserve(per_col);
+  {
+    int i = 0;
+    int j = 0;
+    for(const auto & A : L)
+    {
+      for(int k = 0;k<A.outerSize();++k)
+      {
+        for(typename Eigen::SparseMatrix<Scalar>::InnerIterator it(A,k);it;++it)
+        {
+           Y.insert(i+it.row(),j+k) = it.value();
+        }
+      }
+      i += A.rows();
+      j += A.cols();
+    }
+  }
+}
+
+template <typename DerivedY>
+IGL_INLINE void igl::blkdiag(
+  const std::vector<DerivedY> & L, 
+  Eigen::PlainObjectBase<DerivedY> & Y)
+{
+  int nr = 0;
+  int nc = 0;
+  for(const auto & A : L)
+  {
+    nr += A.rows();
+    nc += A.cols();
+  }
+  Y.setZero(nr,nc);
+  {
+    int i = 0;
+    int j = 0;
+    for(const auto & A : L)
+    {
+      Y.block(i,j,A.rows(),A.cols()) = A;
+      i += A.rows();
+      j += A.cols();
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// explicit template instantiations
+template void igl::blkdiag<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::__1::vector<Eigen::Matrix<double, -1, -1, 0, -1, -1>, std::__1::allocator<Eigen::Matrix<double, -1, -1, 0, -1, -1> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::blkdiag<double>(std::__1::vector<Eigen::SparseMatrix<double, 0, int>, std::__1::allocator<Eigen::SparseMatrix<double, 0, int> > > const&, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/blkdiag.h

@@ -0,0 +1,40 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2020 Alec Jacobson <[email protected]>
+// 
+// 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/.
+#ifndef IGL_BLKDIAG_H
+#define IGL_BLKDIAG_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+#include <vector>
+
+namespace igl
+{
+  // Given a list of matrices place them along the diagonal as blocks of the
+  // output matrix. Like matlab's blkdiag.
+  //
+  // Inputs:
+  //   L  list of matrices {A,B, ...}
+  // Outputs:
+  //   Y  A.rows()+B.rows()+... by A.cols()+B.cols()+... block diagonal
+  //
+  // See also: cat, repdiag
+  template <typename Scalar>
+  IGL_INLINE void blkdiag(
+    const std::vector<Eigen::SparseMatrix<Scalar>> & L, 
+    Eigen::SparseMatrix<Scalar> & Y);
+  template <typename DerivedY>
+  IGL_INLINE void blkdiag(
+    const std::vector<DerivedY> & L, 
+    Eigen::PlainObjectBase<DerivedY> & Y);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "blkdiag.cpp"
+#endif
+
+#endif

tests/include/igl/blkdiag.cpp

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2020 Alec Jacobson <[email protected]>
+// 
+// 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 <test_common.h>
+#include <igl/blkdiag.h>
+
+TEST_CASE("blkdiag: 3-matrices", "[igl]")
+{
+  const Eigen::MatrixXd Ygt= 
+    (Eigen::MatrixXd(3+0+2,3+2+5)<<
+      1,2,3,0,0,0,0,0,0,0,
+      4,5,6,0,0,0,0,0,0,0,
+      7,8,9,0,0,0,0,0,0,0,
+      0,0,0,0,0,0,1,2,3,4,
+      0,0,0,0,0,5,6,7,8,9).finished();
+  Eigen::MatrixXd A(3,3);
+  A<<1,2,3,4,5,6,7,8,9;
+  // This is the correct behavior for a 0×n matrix
+  Eigen::MatrixXd B(0,2);
+  Eigen::MatrixXd C(2,5);
+  C<<0,1,2,3,4,5,6,7,8,9;
+  Eigen::MatrixXd Y;
+  igl::blkdiag({A,B,C},Y);
+  test_common::assert_eq(Y,Ygt);
+  {
+    Eigen::SparseMatrix<double> sA = A.sparseView();
+    Eigen::SparseMatrix<double> sB = B.sparseView();
+    Eigen::SparseMatrix<double> sC = C.sparseView();
+    Eigen::SparseMatrix<double> sY;
+    igl::blkdiag({sA,sB,sC},sY);
+    Y = sY;
+    test_common::assert_eq(Y,Ygt);
+  }
+}
+