1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2011 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 #include "sparse.h" 11 12 template<int OtherStorage, typename SparseMatrixType> void sparse_permutations(const SparseMatrixType& ref) 13 { 14 typedef typename SparseMatrixType::Index Index; 15 16 const Index rows = ref.rows(); 17 const Index cols = ref.cols(); 18 typedef typename SparseMatrixType::Scalar Scalar; 19 typedef typename SparseMatrixType::Index Index; 20 typedef SparseMatrix<Scalar, OtherStorage, Index> OtherSparseMatrixType; 21 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 22 typedef Matrix<Index,Dynamic,1> VectorI; 23 24 double density = (std::max)(8./(rows*cols), 0.01); 25 26 SparseMatrixType mat(rows, cols), up(rows,cols), lo(rows,cols); 27 OtherSparseMatrixType res; 28 DenseMatrix mat_d = DenseMatrix::Zero(rows, cols), up_sym_d, lo_sym_d, res_d; 29 30 initSparse<Scalar>(density, mat_d, mat, 0); 31 32 up = mat.template triangularView<Upper>(); 33 lo = mat.template triangularView<Lower>(); 34 35 up_sym_d = mat_d.template selfadjointView<Upper>(); 36 lo_sym_d = mat_d.template selfadjointView<Lower>(); 37 38 VERIFY_IS_APPROX(mat, mat_d); 39 VERIFY_IS_APPROX(up, DenseMatrix(mat_d.template triangularView<Upper>())); 40 VERIFY_IS_APPROX(lo, DenseMatrix(mat_d.template triangularView<Lower>())); 41 42 PermutationMatrix<Dynamic> p, p_null; 43 VectorI pi; 44 randomPermutationVector(pi, cols); 45 p.indices() = pi; 46 47 res = mat*p; 48 res_d = mat_d*p; 49 VERIFY(res.isApprox(res_d) && "mat*p"); 50 51 res = p*mat; 52 res_d = p*mat_d; 53 VERIFY(res.isApprox(res_d) && "p*mat"); 54 55 res = mat*p.inverse(); 56 res_d = mat*p.inverse(); 57 VERIFY(res.isApprox(res_d) && "mat*inv(p)"); 58 59 res = p.inverse()*mat; 60 res_d = p.inverse()*mat_d; 61 VERIFY(res.isApprox(res_d) && "inv(p)*mat"); 62 63 res = mat.twistedBy(p); 64 res_d = (p * mat_d) * p.inverse(); 65 VERIFY(res.isApprox(res_d) && "p*mat*inv(p)"); 66 67 68 res = mat.template selfadjointView<Upper>().twistedBy(p_null); 69 res_d = up_sym_d; 70 VERIFY(res.isApprox(res_d) && "full selfadjoint upper to full"); 71 72 res = mat.template selfadjointView<Lower>().twistedBy(p_null); 73 res_d = lo_sym_d; 74 VERIFY(res.isApprox(res_d) && "full selfadjoint lower to full"); 75 76 77 res = up.template selfadjointView<Upper>().twistedBy(p_null); 78 res_d = up_sym_d; 79 VERIFY(res.isApprox(res_d) && "upper selfadjoint to full"); 80 81 res = lo.template selfadjointView<Lower>().twistedBy(p_null); 82 res_d = lo_sym_d; 83 VERIFY(res.isApprox(res_d) && "lower selfadjoint full"); 84 85 86 res = mat.template selfadjointView<Upper>(); 87 res_d = up_sym_d; 88 VERIFY(res.isApprox(res_d) && "full selfadjoint upper to full"); 89 90 res = mat.template selfadjointView<Lower>(); 91 res_d = lo_sym_d; 92 VERIFY(res.isApprox(res_d) && "full selfadjoint lower to full"); 93 94 res = up.template selfadjointView<Upper>(); 95 res_d = up_sym_d; 96 VERIFY(res.isApprox(res_d) && "upper selfadjoint to full"); 97 98 res = lo.template selfadjointView<Lower>(); 99 res_d = lo_sym_d; 100 VERIFY(res.isApprox(res_d) && "lower selfadjoint full"); 101 102 103 res.template selfadjointView<Upper>() = mat.template selfadjointView<Upper>(); 104 res_d = up_sym_d.template triangularView<Upper>(); 105 VERIFY(res.isApprox(res_d) && "full selfadjoint upper to upper"); 106 107 res.template selfadjointView<Lower>() = mat.template selfadjointView<Upper>(); 108 res_d = up_sym_d.template triangularView<Lower>(); 109 VERIFY(res.isApprox(res_d) && "full selfadjoint upper to lower"); 110 111 res.template selfadjointView<Upper>() = mat.template selfadjointView<Lower>(); 112 res_d = lo_sym_d.template triangularView<Upper>(); 113 VERIFY(res.isApprox(res_d) && "full selfadjoint lower to upper"); 114 115 res.template selfadjointView<Lower>() = mat.template selfadjointView<Lower>(); 116 res_d = lo_sym_d.template triangularView<Lower>(); 117 VERIFY(res.isApprox(res_d) && "full selfadjoint lower to lower"); 118 119 120 121 res.template selfadjointView<Upper>() = mat.template selfadjointView<Upper>().twistedBy(p); 122 res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView<Upper>(); 123 VERIFY(res.isApprox(res_d) && "full selfadjoint upper twisted to upper"); 124 125 res.template selfadjointView<Upper>() = mat.template selfadjointView<Lower>().twistedBy(p); 126 res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView<Upper>(); 127 VERIFY(res.isApprox(res_d) && "full selfadjoint lower twisted to upper"); 128 129 res.template selfadjointView<Lower>() = mat.template selfadjointView<Lower>().twistedBy(p); 130 res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView<Lower>(); 131 VERIFY(res.isApprox(res_d) && "full selfadjoint lower twisted to lower"); 132 133 res.template selfadjointView<Lower>() = mat.template selfadjointView<Upper>().twistedBy(p); 134 res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView<Lower>(); 135 VERIFY(res.isApprox(res_d) && "full selfadjoint upper twisted to lower"); 136 137 138 res.template selfadjointView<Upper>() = up.template selfadjointView<Upper>().twistedBy(p); 139 res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView<Upper>(); 140 VERIFY(res.isApprox(res_d) && "upper selfadjoint twisted to upper"); 141 142 res.template selfadjointView<Upper>() = lo.template selfadjointView<Lower>().twistedBy(p); 143 res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView<Upper>(); 144 VERIFY(res.isApprox(res_d) && "lower selfadjoint twisted to upper"); 145 146 res.template selfadjointView<Lower>() = lo.template selfadjointView<Lower>().twistedBy(p); 147 res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView<Lower>(); 148 VERIFY(res.isApprox(res_d) && "lower selfadjoint twisted to lower"); 149 150 res.template selfadjointView<Lower>() = up.template selfadjointView<Upper>().twistedBy(p); 151 res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView<Lower>(); 152 VERIFY(res.isApprox(res_d) && "upper selfadjoint twisted to lower"); 153 154 155 res = mat.template selfadjointView<Upper>().twistedBy(p); 156 res_d = (p * up_sym_d) * p.inverse(); 157 VERIFY(res.isApprox(res_d) && "full selfadjoint upper twisted to full"); 158 159 res = mat.template selfadjointView<Lower>().twistedBy(p); 160 res_d = (p * lo_sym_d) * p.inverse(); 161 VERIFY(res.isApprox(res_d) && "full selfadjoint lower twisted to full"); 162 163 res = up.template selfadjointView<Upper>().twistedBy(p); 164 res_d = (p * up_sym_d) * p.inverse(); 165 VERIFY(res.isApprox(res_d) && "upper selfadjoint twisted to full"); 166 167 res = lo.template selfadjointView<Lower>().twistedBy(p); 168 res_d = (p * lo_sym_d) * p.inverse(); 169 VERIFY(res.isApprox(res_d) && "lower selfadjoint twisted to full"); 170 } 171 172 template<typename Scalar> void sparse_permutations_all(int size) 173 { 174 CALL_SUBTEST(( sparse_permutations<ColMajor>(SparseMatrix<Scalar, ColMajor>(size,size)) )); 175 CALL_SUBTEST(( sparse_permutations<ColMajor>(SparseMatrix<Scalar, RowMajor>(size,size)) )); 176 CALL_SUBTEST(( sparse_permutations<RowMajor>(SparseMatrix<Scalar, ColMajor>(size,size)) )); 177 CALL_SUBTEST(( sparse_permutations<RowMajor>(SparseMatrix<Scalar, RowMajor>(size,size)) )); 178 } 179 180 void test_sparse_permutations() 181 { 182 for(int i = 0; i < g_repeat; i++) { 183 int s = Eigen::internal::random<int>(1,50); 184 CALL_SUBTEST_1(( sparse_permutations_all<double>(s) )); 185 CALL_SUBTEST_2(( sparse_permutations_all<std::complex<double> >(s) )); 186 } 187 } 188
