ground/gcs/src/libs/eigen/test/sparse_permutations.cpp

   1 // This file is part of Eigen, a lightweight C++ template library
   2 // for linear algebra.
   3 //
   4 // Copyright (C) 2011-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
   5 //
   6 // This Source Code Form is subject to the terms of the Mozilla
   7 // Public License v. 2.0. If a copy of the MPL was not distributed
   8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
   9
  10
  11 static long int nb_transposed_copies;
  12 #define EIGEN_SPARSE_TRANSPOSED_COPY_PLUGIN {nb_transposed_copies++;}
  13 #define VERIFY_TRANSPOSITION_COUNT(XPR,N) {\
  14     nb_transposed_copies = 0; \
  15     XPR; \
  16     if(nb_transposed_copies!=N) std::cerr << "nb_transposed_copies == " << nb_transposed_copies << "\n"; \
  17     VERIFY( (#XPR) && nb_transposed_copies==N ); \
  18   }
  19
  20 #include "sparse.h"
  21
  22 template<typename T>
  23 bool is_sorted(const T& mat) {
  24   for(Index k = 0; k<mat.outerSize(); ++k)
  25   {
  26     Index prev = -1;
  27     for(typename T::InnerIterator it(mat,k); it; ++it)
  28     {
  29       if(prev>=it.index())
  30         return false;
  31       prev = it.index();
  32     }
  33   }
  34   return true;
  35 }
  36
  37 template<typename T>
  38 typename internal::nested_eval<T,1>::type eval(const T &xpr)
  39 {
  40   VERIFY( int(internal::nested_eval<T,1>::type::Flags&RowMajorBit) == int(internal::evaluator<T>::Flags&RowMajorBit) );
  41   return xpr;
  42 }
  43
  44 template<int OtherStorage, typename SparseMatrixType> void sparse_permutations(const SparseMatrixType& ref)
  45 {
  46   const Index rows = ref.rows();
  47   const Index cols = ref.cols();
  48   typedef typename SparseMatrixType::Scalar Scalar;
  49   typedef typename SparseMatrixType::StorageIndex StorageIndex;
  50   typedef SparseMatrix<Scalar, OtherStorage, StorageIndex> OtherSparseMatrixType;
  51   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
  52   typedef Matrix<StorageIndex,Dynamic,1> VectorI;
  53 //   bool IsRowMajor1 = SparseMatrixType::IsRowMajor;
  54 //   bool IsRowMajor2 = OtherSparseMatrixType::IsRowMajor;
  55
  56   double density = (std::max)(8./(rows*cols), 0.01);
  57
  58   SparseMatrixType mat(rows, cols), up(rows,cols), lo(rows,cols);
  59   OtherSparseMatrixType res;
  60   DenseMatrix mat_d = DenseMatrix::Zero(rows, cols), up_sym_d, lo_sym_d, res_d;
  61
  62   initSparse<Scalar>(density, mat_d, mat, 0);
  63
  64   up = mat.template triangularView<Upper>();
  65   lo = mat.template triangularView<Lower>();
  66
  67   up_sym_d = mat_d.template selfadjointView<Upper>();
  68   lo_sym_d = mat_d.template selfadjointView<Lower>();
  69
  70   VERIFY_IS_APPROX(mat, mat_d);
  71   VERIFY_IS_APPROX(up, DenseMatrix(mat_d.template triangularView<Upper>()));
  72   VERIFY_IS_APPROX(lo, DenseMatrix(mat_d.template triangularView<Lower>()));
  73
  74   PermutationMatrix<Dynamic> p, p_null;
  75   VectorI pi;
  76   randomPermutationVector(pi, cols);
  77   p.indices() = pi;
  78
  79   VERIFY( is_sorted( ::eval(mat*p) ));
  80   VERIFY( is_sorted( res = mat*p ));
  81   VERIFY_TRANSPOSITION_COUNT( ::eval(mat*p), 0);
  82   //VERIFY_TRANSPOSITION_COUNT( res = mat*p, IsRowMajor ? 1 : 0 );
  83   res_d = mat_d*p;
  84   VERIFY(res.isApprox(res_d) && "mat*p");
  85
  86   VERIFY( is_sorted( ::eval(p*mat) ));
  87   VERIFY( is_sorted( res = p*mat ));
  88   VERIFY_TRANSPOSITION_COUNT( ::eval(p*mat), 0);
  89   res_d = p*mat_d;
  90   VERIFY(res.isApprox(res_d) && "p*mat");
  91
  92   VERIFY( is_sorted( (mat*p).eval() ));
  93   VERIFY( is_sorted( res = mat*p.inverse() ));
  94   VERIFY_TRANSPOSITION_COUNT( ::eval(mat*p.inverse()), 0);
  95   res_d = mat*p.inverse();
  96   VERIFY(res.isApprox(res_d) && "mat*inv(p)");
  97
  98   VERIFY( is_sorted( (p*mat+p*mat).eval() ));
  99   VERIFY( is_sorted( res = p.inverse()*mat ));
 100   VERIFY_TRANSPOSITION_COUNT( ::eval(p.inverse()*mat), 0);
 101   res_d = p.inverse()*mat_d;
 102   VERIFY(res.isApprox(res_d) && "inv(p)*mat");
 103
 104   VERIFY( is_sorted( (p * mat * p.inverse()).eval() ));
 105   VERIFY( is_sorted( res = mat.twistedBy(p) ));
 106   VERIFY_TRANSPOSITION_COUNT( ::eval(p * mat * p.inverse()), 0);
 107   res_d = (p * mat_d) * p.inverse();
 108   VERIFY(res.isApprox(res_d) && "p*mat*inv(p)");
 109
 110
 111   VERIFY( is_sorted( res = mat.template selfadjointView<Upper>().twistedBy(p_null) ));
 112   res_d = up_sym_d;
 113   VERIFY(res.isApprox(res_d) && "full selfadjoint upper to full");
 114
 115   VERIFY( is_sorted( res = mat.template selfadjointView<Lower>().twistedBy(p_null) ));
 116   res_d = lo_sym_d;
 117   VERIFY(res.isApprox(res_d) && "full selfadjoint lower to full");
 118
 119
 120   VERIFY( is_sorted( res = up.template selfadjointView<Upper>().twistedBy(p_null) ));
 121   res_d = up_sym_d;
 122   VERIFY(res.isApprox(res_d) && "upper selfadjoint to full");
 123
 124   VERIFY( is_sorted( res = lo.template selfadjointView<Lower>().twistedBy(p_null) ));
 125   res_d = lo_sym_d;
 126   VERIFY(res.isApprox(res_d) && "lower selfadjoint full");
 127
 128
 129   VERIFY( is_sorted( res = mat.template selfadjointView<Upper>() ));
 130   res_d = up_sym_d;
 131   VERIFY(res.isApprox(res_d) && "full selfadjoint upper to full");
 132
 133   VERIFY( is_sorted( res = mat.template selfadjointView<Lower>() ));
 134   res_d = lo_sym_d;
 135   VERIFY(res.isApprox(res_d) && "full selfadjoint lower to full");
 136
 137   VERIFY( is_sorted( res = up.template selfadjointView<Upper>() ));
 138   res_d = up_sym_d;
 139   VERIFY(res.isApprox(res_d) && "upper selfadjoint to full");
 140
 141   VERIFY( is_sorted( res = lo.template selfadjointView<Lower>() ));
 142   res_d = lo_sym_d;
 143   VERIFY(res.isApprox(res_d) && "lower selfadjoint full");
 144
 145
 146   res.template selfadjointView<Upper>() = mat.template selfadjointView<Upper>();
 147   res_d = up_sym_d.template triangularView<Upper>();
 148   VERIFY(res.isApprox(res_d) && "full selfadjoint upper to upper");
 149
 150   res.template selfadjointView<Lower>() = mat.template selfadjointView<Upper>();
 151   res_d = up_sym_d.template triangularView<Lower>();
 152   VERIFY(res.isApprox(res_d) && "full selfadjoint upper to lower");
 153
 154   res.template selfadjointView<Upper>() = mat.template selfadjointView<Lower>();
 155   res_d = lo_sym_d.template triangularView<Upper>();
 156   VERIFY(res.isApprox(res_d) && "full selfadjoint lower to upper");
 157
 158   res.template selfadjointView<Lower>() = mat.template selfadjointView<Lower>();
 159   res_d = lo_sym_d.template triangularView<Lower>();
 160   VERIFY(res.isApprox(res_d) && "full selfadjoint lower to lower");
 161
 162
 163
 164   res.template selfadjointView<Upper>() = mat.template selfadjointView<Upper>().twistedBy(p);
 165   res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView<Upper>();
 166   VERIFY(res.isApprox(res_d) && "full selfadjoint upper twisted to upper");
 167
 168   res.template selfadjointView<Upper>() = mat.template selfadjointView<Lower>().twistedBy(p);
 169   res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView<Upper>();
 170   VERIFY(res.isApprox(res_d) && "full selfadjoint lower twisted to upper");
 171
 172   res.template selfadjointView<Lower>() = mat.template selfadjointView<Lower>().twistedBy(p);
 173   res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView<Lower>();
 174   VERIFY(res.isApprox(res_d) && "full selfadjoint lower twisted to lower");
 175
 176   res.template selfadjointView<Lower>() = mat.template selfadjointView<Upper>().twistedBy(p);
 177   res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView<Lower>();
 178   VERIFY(res.isApprox(res_d) && "full selfadjoint upper twisted to lower");
 179
 180
 181   res.template selfadjointView<Upper>() = up.template selfadjointView<Upper>().twistedBy(p);
 182   res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView<Upper>();
 183   VERIFY(res.isApprox(res_d) && "upper selfadjoint twisted to upper");
 184
 185   res.template selfadjointView<Upper>() = lo.template selfadjointView<Lower>().twistedBy(p);
 186   res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView<Upper>();
 187   VERIFY(res.isApprox(res_d) && "lower selfadjoint twisted to upper");
 188
 189   res.template selfadjointView<Lower>() = lo.template selfadjointView<Lower>().twistedBy(p);
 190   res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView<Lower>();
 191   VERIFY(res.isApprox(res_d) && "lower selfadjoint twisted to lower");
 192
 193   res.template selfadjointView<Lower>() = up.template selfadjointView<Upper>().twistedBy(p);
 194   res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView<Lower>();
 195   VERIFY(res.isApprox(res_d) && "upper selfadjoint twisted to lower");
 196
 197
 198   VERIFY( is_sorted( res = mat.template selfadjointView<Upper>().twistedBy(p) ));
 199   res_d = (p * up_sym_d) * p.inverse();
 200   VERIFY(res.isApprox(res_d) && "full selfadjoint upper twisted to full");
 201
 202   VERIFY( is_sorted( res = mat.template selfadjointView<Lower>().twistedBy(p) ));
 203   res_d = (p * lo_sym_d) * p.inverse();
 204   VERIFY(res.isApprox(res_d) && "full selfadjoint lower twisted to full");
 205
 206   VERIFY( is_sorted( res = up.template selfadjointView<Upper>().twistedBy(p) ));
 207   res_d = (p * up_sym_d) * p.inverse();
 208   VERIFY(res.isApprox(res_d) && "upper selfadjoint twisted to full");
 209
 210   VERIFY( is_sorted( res = lo.template selfadjointView<Lower>().twistedBy(p) ));
 211   res_d = (p * lo_sym_d) * p.inverse();
 212   VERIFY(res.isApprox(res_d) && "lower selfadjoint twisted to full");
 213 }
 214
 215 template<typename Scalar> void sparse_permutations_all(int size)
 216 {
 217   CALL_SUBTEST(( sparse_permutations<ColMajor>(SparseMatrix<Scalar, ColMajor>(size,size)) ));
 218   CALL_SUBTEST(( sparse_permutations<ColMajor>(SparseMatrix<Scalar, RowMajor>(size,size)) ));
 219   CALL_SUBTEST(( sparse_permutations<RowMajor>(SparseMatrix<Scalar, ColMajor>(size,size)) ));
 220   CALL_SUBTEST(( sparse_permutations<RowMajor>(SparseMatrix<Scalar, RowMajor>(size,size)) ));
 221 }
 222
 223 void test_sparse_permutations()
 224 {
 225   for(int i = 0; i < g_repeat; i++) {
 226     int s = Eigen::internal::random<int>(1,50);
 227     CALL_SUBTEST_1((  sparse_permutations_all<double>(s) ));
 228     CALL_SUBTEST_2((  sparse_permutations_all<std::complex<double> >(s) ));
 229   }
 230
 231   VERIFY((internal::is_same<internal::permutation_matrix_product<SparseMatrix<double>,OnTheRight,false,SparseShape>::ReturnType,
 232                             internal::nested_eval<Product<SparseMatrix<double>,PermutationMatrix<Dynamic,Dynamic>,AliasFreeProduct>,1>::type>::value));
 233
 234   VERIFY((internal::is_same<internal::permutation_matrix_product<SparseMatrix<double>,OnTheLeft,false,SparseShape>::ReturnType,
 235                             internal::nested_eval<Product<PermutationMatrix<Dynamic,Dynamic>,SparseMatrix<double>,AliasFreeProduct>,1>::type>::value));
 236 }