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 (at) 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 } 237
