Skip to content
Normal view History Permalink
sparse_permutations.cpp 9.75 KiB
Newer Older
Luker's avatar
CMake rewrite
Luker committed
1 2 3 
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
4 
// Copyright (C) 2011-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
Luker's avatar
CMake rewrite
Luker committed
5 6 7 8 9 
//
// 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/.
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
10 11 12 13 14 15 16 17 18 19 

static long int nb_transposed_copies;
#define EIGEN_SPARSE_TRANSPOSED_COPY_PLUGIN {nb_transposed_copies++;}
#define VERIFY_TRANSPOSITION_COUNT(XPR,N) {\
    nb_transposed_copies = 0; \
    XPR; \
    if(nb_transposed_copies!=N) std::cerr << "nb_transposed_copies == " << nb_transposed_copies << "\n"; \
    VERIFY( (#XPR) && nb_transposed_copies==N ); \
  }
Luker's avatar
CMake rewrite
Luker committed
20 21 
#include "sparse.h"
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 
template<typename T>
bool is_sorted(const T& mat) {
  for(Index k = 0; k<mat.outerSize(); ++k)
  {
    Index prev = -1;
    for(typename T::InnerIterator it(mat,k); it; ++it)
    {
      if(prev>=it.index())
        return false;
      prev = it.index();
    }
  }
  return true;
}

template<typename T>
typename internal::nested_eval<T,1>::type eval(const T &xpr)
Luker's avatar
CMake rewrite
Luker committed
39 
{
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
40 41 42 
  VERIFY( int(internal::nested_eval<T,1>::type::Flags&RowMajorBit) == int(internal::evaluator<T>::Flags&RowMajorBit) );
  return xpr;
}
Luker's avatar
CMake rewrite
Luker committed
43
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
44 45 
template<int OtherStorage, typename SparseMatrixType> void sparse_permutations(const SparseMatrixType& ref)
{
Luker's avatar
CMake rewrite
Luker committed
46 47 48 
  const Index rows = ref.rows();
  const Index cols = ref.cols();
  typedef typename SparseMatrixType::Scalar Scalar;
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
49 50 
  typedef typename SparseMatrixType::StorageIndex StorageIndex;
  typedef SparseMatrix<Scalar, OtherStorage, StorageIndex> OtherSparseMatrixType;
Luker's avatar
CMake rewrite
Luker committed
51 
  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
52 53 54 
  typedef Matrix<StorageIndex,Dynamic,1> VectorI;
//   bool IsRowMajor1 = SparseMatrixType::IsRowMajor;
//   bool IsRowMajor2 = OtherSparseMatrixType::IsRowMajor;
Luker's avatar
CMake rewrite
Luker committed
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 
  
  double density = (std::max)(8./(rows*cols), 0.01);
  
  SparseMatrixType mat(rows, cols), up(rows,cols), lo(rows,cols);
  OtherSparseMatrixType res;
  DenseMatrix mat_d = DenseMatrix::Zero(rows, cols), up_sym_d, lo_sym_d, res_d;
  
  initSparse<Scalar>(density, mat_d, mat, 0);

  up = mat.template triangularView<Upper>();
  lo = mat.template triangularView<Lower>();
  
  up_sym_d = mat_d.template selfadjointView<Upper>();
  lo_sym_d = mat_d.template selfadjointView<Lower>();
  
  VERIFY_IS_APPROX(mat, mat_d);
  VERIFY_IS_APPROX(up, DenseMatrix(mat_d.template triangularView<Upper>()));
  VERIFY_IS_APPROX(lo, DenseMatrix(mat_d.template triangularView<Lower>()));
  
  PermutationMatrix<Dynamic> p, p_null;
  VectorI pi;
  randomPermutationVector(pi, cols);
  p.indices() = pi;
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
79 80 81 82 
  VERIFY( is_sorted( ::eval(mat*p) ));
  VERIFY( is_sorted( res = mat*p ));
  VERIFY_TRANSPOSITION_COUNT( ::eval(mat*p), 0);
  //VERIFY_TRANSPOSITION_COUNT( res = mat*p, IsRowMajor ? 1 : 0 );
Luker's avatar
CMake rewrite
Luker committed
83 84 85 
  res_d = mat_d*p;
  VERIFY(res.isApprox(res_d) && "mat*p");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
86 87 88 
  VERIFY( is_sorted( ::eval(p*mat) ));
  VERIFY( is_sorted( res = p*mat ));
  VERIFY_TRANSPOSITION_COUNT( ::eval(p*mat), 0);
Luker's avatar
CMake rewrite
Luker committed
89 90 91 
  res_d = p*mat_d;
  VERIFY(res.isApprox(res_d) && "p*mat");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
92 93 94 
  VERIFY( is_sorted( (mat*p).eval() ));
  VERIFY( is_sorted( res = mat*p.inverse() ));
  VERIFY_TRANSPOSITION_COUNT( ::eval(mat*p.inverse()), 0);
Luker's avatar
CMake rewrite
Luker committed
95 96 97 
  res_d = mat*p.inverse();
  VERIFY(res.isApprox(res_d) && "mat*inv(p)");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
98 99 100 
  VERIFY( is_sorted( (p*mat+p*mat).eval() ));
  VERIFY( is_sorted( res = p.inverse()*mat ));
  VERIFY_TRANSPOSITION_COUNT( ::eval(p.inverse()*mat), 0);
Luker's avatar
CMake rewrite
Luker committed
101 102 103 
  res_d = p.inverse()*mat_d;
  VERIFY(res.isApprox(res_d) && "inv(p)*mat");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
104 105 106 
  VERIFY( is_sorted( (p * mat * p.inverse()).eval() ));
  VERIFY( is_sorted( res = mat.twistedBy(p) ));
  VERIFY_TRANSPOSITION_COUNT( ::eval(p * mat * p.inverse()), 0);
Luker's avatar
CMake rewrite
Luker committed
107 108 109 110 
  res_d = (p * mat_d) * p.inverse();
  VERIFY(res.isApprox(res_d) && "p*mat*inv(p)");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
111 
  VERIFY( is_sorted( res = mat.template selfadjointView<Upper>().twistedBy(p_null) ));
Luker's avatar
CMake rewrite
Luker committed
112 113 114 
  res_d = up_sym_d;
  VERIFY(res.isApprox(res_d) && "full selfadjoint upper to full");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
115 
  VERIFY( is_sorted( res = mat.template selfadjointView<Lower>().twistedBy(p_null) ));
Luker's avatar
CMake rewrite
Luker committed
116 117 118 119 
  res_d = lo_sym_d;
  VERIFY(res.isApprox(res_d) && "full selfadjoint lower to full");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
120 
  VERIFY( is_sorted( res = up.template selfadjointView<Upper>().twistedBy(p_null) ));
Luker's avatar
CMake rewrite
Luker committed
121 122 123 
  res_d = up_sym_d;
  VERIFY(res.isApprox(res_d) && "upper selfadjoint to full");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
124 
  VERIFY( is_sorted( res = lo.template selfadjointView<Lower>().twistedBy(p_null) ));
Luker's avatar
CMake rewrite
Luker committed
125 126 127 128 
  res_d = lo_sym_d;
  VERIFY(res.isApprox(res_d) && "lower selfadjoint full");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
129 
  VERIFY( is_sorted( res = mat.template selfadjointView<Upper>() ));
Luker's avatar
CMake rewrite
Luker committed
130 131 132 
  res_d = up_sym_d;
  VERIFY(res.isApprox(res_d) && "full selfadjoint upper to full");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
133 
  VERIFY( is_sorted( res = mat.template selfadjointView<Lower>() ));
Luker's avatar
CMake rewrite
Luker committed
134 135 136 
  res_d = lo_sym_d;
  VERIFY(res.isApprox(res_d) && "full selfadjoint lower to full");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
137 
  VERIFY( is_sorted( res = up.template selfadjointView<Upper>() ));
Luker's avatar
CMake rewrite
Luker committed
138 139 140 
  res_d = up_sym_d;
  VERIFY(res.isApprox(res_d) && "upper selfadjoint to full");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
141 
  VERIFY( is_sorted( res = lo.template selfadjointView<Lower>() ));
Luker's avatar
CMake rewrite
Luker committed
142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 
  res_d = lo_sym_d;
  VERIFY(res.isApprox(res_d) && "lower selfadjoint full");


  res.template selfadjointView<Upper>() = mat.template selfadjointView<Upper>();
  res_d = up_sym_d.template triangularView<Upper>();
  VERIFY(res.isApprox(res_d) && "full selfadjoint upper to upper");

  res.template selfadjointView<Lower>() = mat.template selfadjointView<Upper>();
  res_d = up_sym_d.template triangularView<Lower>();
  VERIFY(res.isApprox(res_d) && "full selfadjoint upper to lower");

  res.template selfadjointView<Upper>() = mat.template selfadjointView<Lower>();
  res_d = lo_sym_d.template triangularView<Upper>();
  VERIFY(res.isApprox(res_d) && "full selfadjoint lower to upper");

  res.template selfadjointView<Lower>() = mat.template selfadjointView<Lower>();
  res_d = lo_sym_d.template triangularView<Lower>();
  VERIFY(res.isApprox(res_d) && "full selfadjoint lower to lower");

  
  
  res.template selfadjointView<Upper>() = mat.template selfadjointView<Upper>().twistedBy(p);
  res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView<Upper>();
  VERIFY(res.isApprox(res_d) && "full selfadjoint upper twisted to upper");
  
  res.template selfadjointView<Upper>() = mat.template selfadjointView<Lower>().twistedBy(p);
  res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView<Upper>();
  VERIFY(res.isApprox(res_d) && "full selfadjoint lower twisted to upper");
  
  res.template selfadjointView<Lower>() = mat.template selfadjointView<Lower>().twistedBy(p);
  res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView<Lower>();
  VERIFY(res.isApprox(res_d) && "full selfadjoint lower twisted to lower");
  
  res.template selfadjointView<Lower>() = mat.template selfadjointView<Upper>().twistedBy(p);
  res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView<Lower>();
  VERIFY(res.isApprox(res_d) && "full selfadjoint upper twisted to lower");
  
  
  res.template selfadjointView<Upper>() = up.template selfadjointView<Upper>().twistedBy(p);
  res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView<Upper>();
  VERIFY(res.isApprox(res_d) && "upper selfadjoint twisted to upper");
  
  res.template selfadjointView<Upper>() = lo.template selfadjointView<Lower>().twistedBy(p);
  res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView<Upper>();
  VERIFY(res.isApprox(res_d) && "lower selfadjoint twisted to upper");
  
  res.template selfadjointView<Lower>() = lo.template selfadjointView<Lower>().twistedBy(p);
  res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView<Lower>();
  VERIFY(res.isApprox(res_d) && "lower selfadjoint twisted to lower");
  
  res.template selfadjointView<Lower>() = up.template selfadjointView<Upper>().twistedBy(p);
  res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView<Lower>();
  VERIFY(res.isApprox(res_d) && "upper selfadjoint twisted to lower");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
198 
  VERIFY( is_sorted( res = mat.template selfadjointView<Upper>().twistedBy(p) ));
Luker's avatar
CMake rewrite
Luker committed
199 200 201 
  res_d = (p * up_sym_d) * p.inverse();
  VERIFY(res.isApprox(res_d) && "full selfadjoint upper twisted to full");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
202 
  VERIFY( is_sorted( res = mat.template selfadjointView<Lower>().twistedBy(p) ));
Luker's avatar
CMake rewrite
Luker committed
203 204 205 
  res_d = (p * lo_sym_d) * p.inverse();
  VERIFY(res.isApprox(res_d) && "full selfadjoint lower twisted to full");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
206 
  VERIFY( is_sorted( res = up.template selfadjointView<Upper>().twistedBy(p) ));
Luker's avatar
CMake rewrite
Luker committed
207 208 209 
  res_d = (p * up_sym_d) * p.inverse();
  VERIFY(res.isApprox(res_d) && "upper selfadjoint twisted to full");
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
210 
  VERIFY( is_sorted( res = lo.template selfadjointView<Lower>().twistedBy(p) ));
Luker's avatar
CMake rewrite
Luker committed
211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 
  res_d = (p * lo_sym_d) * p.inverse();
  VERIFY(res.isApprox(res_d) && "lower selfadjoint twisted to full");
}

template<typename Scalar> void sparse_permutations_all(int size)
{
  CALL_SUBTEST(( sparse_permutations<ColMajor>(SparseMatrix<Scalar, ColMajor>(size,size)) ));
  CALL_SUBTEST(( sparse_permutations<ColMajor>(SparseMatrix<Scalar, RowMajor>(size,size)) ));
  CALL_SUBTEST(( sparse_permutations<RowMajor>(SparseMatrix<Scalar, ColMajor>(size,size)) ));
  CALL_SUBTEST(( sparse_permutations<RowMajor>(SparseMatrix<Scalar, RowMajor>(size,size)) ));
}

void test_sparse_permutations()
{
  for(int i = 0; i < g_repeat; i++) {
    int s = Eigen::internal::random<int>(1,50);
    CALL_SUBTEST_1((  sparse_permutations_all<double>(s) ));
    CALL_SUBTEST_2((  sparse_permutations_all<std::complex<double> >(s) ));
  }
Luker's avatar
Update to Eigen 3.3.4  
Luker committed
230 231 232 233 234 235 

  VERIFY((internal::is_same<internal::permutation_matrix_product<SparseMatrix<double>,OnTheRight,false,SparseShape>::ReturnType,
                            internal::nested_eval<Product<SparseMatrix<double>,PermutationMatrix<Dynamic,Dynamic>,AliasFreeProduct>,1>::type>::value));

  VERIFY((internal::is_same<internal::permutation_matrix_product<SparseMatrix<double>,OnTheLeft,false,SparseShape>::ReturnType,
                            internal::nested_eval<Product<PermutationMatrix<Dynamic,Dynamic>,SparseMatrix<double>,AliasFreeProduct>,1>::type>::value));
Luker's avatar
CMake rewrite
Luker committed
236 
}