triangular.cpp (11918B)
1 // This file is triangularView of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2009 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 #ifdef EIGEN_TEST_PART_100 11 # define EIGEN_NO_DEPRECATED_WARNING 12 #endif 13 14 #include "main.h" 15 16 17 template<typename MatrixType> void triangular_deprecated(const MatrixType &m) 18 { 19 Index rows = m.rows(); 20 Index cols = m.cols(); 21 MatrixType m1, m2, m3, m4; 22 m1.setRandom(rows,cols); 23 m2.setRandom(rows,cols); 24 m3 = m1; m4 = m2; 25 // deprecated method: 26 m1.template triangularView<Eigen::Upper>().swap(m2); 27 // use this method instead: 28 m3.template triangularView<Eigen::Upper>().swap(m4.template triangularView<Eigen::Upper>()); 29 VERIFY_IS_APPROX(m1,m3); 30 VERIFY_IS_APPROX(m2,m4); 31 // deprecated method: 32 m1.template triangularView<Eigen::Lower>().swap(m4); 33 // use this method instead: 34 m3.template triangularView<Eigen::Lower>().swap(m2.template triangularView<Eigen::Lower>()); 35 VERIFY_IS_APPROX(m1,m3); 36 VERIFY_IS_APPROX(m2,m4); 37 } 38 39 40 template<typename MatrixType> void triangular_square(const MatrixType& m) 41 { 42 typedef typename MatrixType::Scalar Scalar; 43 typedef typename NumTraits<Scalar>::Real RealScalar; 44 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 45 46 RealScalar largerEps = 10*test_precision<RealScalar>(); 47 48 Index rows = m.rows(); 49 Index cols = m.cols(); 50 51 MatrixType m1 = MatrixType::Random(rows, cols), 52 m2 = MatrixType::Random(rows, cols), 53 m3(rows, cols), 54 m4(rows, cols), 55 r1(rows, cols), 56 r2(rows, cols); 57 VectorType v2 = VectorType::Random(rows); 58 59 MatrixType m1up = m1.template triangularView<Upper>(); 60 MatrixType m2up = m2.template triangularView<Upper>(); 61 62 if (rows*cols>1) 63 { 64 VERIFY(m1up.isUpperTriangular()); 65 VERIFY(m2up.transpose().isLowerTriangular()); 66 VERIFY(!m2.isLowerTriangular()); 67 } 68 69 // VERIFY_IS_APPROX(m1up.transpose() * m2, m1.upper().transpose().lower() * m2); 70 71 // test overloaded operator+= 72 r1.setZero(); 73 r2.setZero(); 74 r1.template triangularView<Upper>() += m1; 75 r2 += m1up; 76 VERIFY_IS_APPROX(r1,r2); 77 78 // test overloaded operator= 79 m1.setZero(); 80 m1.template triangularView<Upper>() = m2.transpose() + m2; 81 m3 = m2.transpose() + m2; 82 VERIFY_IS_APPROX(m3.template triangularView<Lower>().transpose().toDenseMatrix(), m1); 83 84 // test overloaded operator= 85 m1.setZero(); 86 m1.template triangularView<Lower>() = m2.transpose() + m2; 87 VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1); 88 89 VERIFY_IS_APPROX(m3.template triangularView<Lower>().conjugate().toDenseMatrix(), 90 m3.conjugate().template triangularView<Lower>().toDenseMatrix()); 91 92 m1 = MatrixType::Random(rows, cols); 93 for (int i=0; i<rows; ++i) 94 while (numext::abs2(m1(i,i))<RealScalar(1e-1)) m1(i,i) = internal::random<Scalar>(); 95 96 Transpose<MatrixType> trm4(m4); 97 // test back and forward substitution with a vector as the rhs 98 m3 = m1.template triangularView<Upper>(); 99 VERIFY(v2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(v2)), largerEps)); 100 m3 = m1.template triangularView<Lower>(); 101 VERIFY(v2.isApprox(m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(v2)), largerEps)); 102 m3 = m1.template triangularView<Upper>(); 103 VERIFY(v2.isApprox(m3 * (m1.template triangularView<Upper>().solve(v2)), largerEps)); 104 m3 = m1.template triangularView<Lower>(); 105 VERIFY(v2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(v2)), largerEps)); 106 107 // test back and forward substitution with a matrix as the rhs 108 m3 = m1.template triangularView<Upper>(); 109 VERIFY(m2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(m2)), largerEps)); 110 m3 = m1.template triangularView<Lower>(); 111 VERIFY(m2.isApprox(m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(m2)), largerEps)); 112 m3 = m1.template triangularView<Upper>(); 113 VERIFY(m2.isApprox(m3 * (m1.template triangularView<Upper>().solve(m2)), largerEps)); 114 m3 = m1.template triangularView<Lower>(); 115 VERIFY(m2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(m2)), largerEps)); 116 117 // check M * inv(L) using in place API 118 m4 = m3; 119 m1.transpose().template triangularView<Eigen::Upper>().solveInPlace(trm4); 120 VERIFY_IS_APPROX(m4 * m1.template triangularView<Eigen::Lower>(), m3); 121 122 // check M * inv(U) using in place API 123 m3 = m1.template triangularView<Upper>(); 124 m4 = m3; 125 m3.transpose().template triangularView<Eigen::Lower>().solveInPlace(trm4); 126 VERIFY_IS_APPROX(m4 * m1.template triangularView<Eigen::Upper>(), m3); 127 128 // check solve with unit diagonal 129 m3 = m1.template triangularView<UnitUpper>(); 130 VERIFY(m2.isApprox(m3 * (m1.template triangularView<UnitUpper>().solve(m2)), largerEps)); 131 132 // VERIFY(( m1.template triangularView<Upper>() 133 // * m2.template triangularView<Upper>()).isUpperTriangular()); 134 135 // test swap 136 m1.setOnes(); 137 m2.setZero(); 138 m2.template triangularView<Upper>().swap(m1.template triangularView<Eigen::Upper>()); 139 m3.setZero(); 140 m3.template triangularView<Upper>().setOnes(); 141 VERIFY_IS_APPROX(m2,m3); 142 VERIFY_RAISES_STATIC_ASSERT(m1.template triangularView<Eigen::Lower>().swap(m2.template triangularView<Eigen::Upper>())); 143 144 m1.setRandom(); 145 m3 = m1.template triangularView<Upper>(); 146 Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic> m5(cols, internal::random<int>(1,20)); m5.setRandom(); 147 Matrix<Scalar, Dynamic, MatrixType::RowsAtCompileTime> m6(internal::random<int>(1,20), rows); m6.setRandom(); 148 VERIFY_IS_APPROX(m1.template triangularView<Upper>() * m5, m3*m5); 149 VERIFY_IS_APPROX(m6*m1.template triangularView<Upper>(), m6*m3); 150 151 m1up = m1.template triangularView<Upper>(); 152 VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().template triangularView<Upper>().toDenseMatrix(), m1up); 153 VERIFY_IS_APPROX(m1up.template selfadjointView<Upper>().template triangularView<Upper>().toDenseMatrix(), m1up); 154 VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().template triangularView<Lower>().toDenseMatrix(), m1up.adjoint()); 155 VERIFY_IS_APPROX(m1up.template selfadjointView<Upper>().template triangularView<Lower>().toDenseMatrix(), m1up.adjoint()); 156 157 VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().diagonal(), m1.diagonal()); 158 159 m3.setRandom(); 160 const MatrixType& m3c(m3); 161 VERIFY( is_same_type(m3c.template triangularView<Lower>(),m3.template triangularView<Lower>().template conjugateIf<false>()) ); 162 VERIFY( is_same_type(m3c.template triangularView<Lower>().conjugate(),m3.template triangularView<Lower>().template conjugateIf<true>()) ); 163 VERIFY_IS_APPROX(m3.template triangularView<Lower>().template conjugateIf<true>().toDenseMatrix(), 164 m3.conjugate().template triangularView<Lower>().toDenseMatrix()); 165 VERIFY_IS_APPROX(m3.template triangularView<Lower>().template conjugateIf<false>().toDenseMatrix(), 166 m3.template triangularView<Lower>().toDenseMatrix()); 167 168 VERIFY( is_same_type(m3c.template selfadjointView<Lower>(),m3.template selfadjointView<Lower>().template conjugateIf<false>()) ); 169 VERIFY( is_same_type(m3c.template selfadjointView<Lower>().conjugate(),m3.template selfadjointView<Lower>().template conjugateIf<true>()) ); 170 VERIFY_IS_APPROX(m3.template selfadjointView<Lower>().template conjugateIf<true>().toDenseMatrix(), 171 m3.conjugate().template selfadjointView<Lower>().toDenseMatrix()); 172 VERIFY_IS_APPROX(m3.template selfadjointView<Lower>().template conjugateIf<false>().toDenseMatrix(), 173 m3.template selfadjointView<Lower>().toDenseMatrix()); 174 175 } 176 177 178 template<typename MatrixType> void triangular_rect(const MatrixType& m) 179 { 180 typedef typename MatrixType::Scalar Scalar; 181 typedef typename NumTraits<Scalar>::Real RealScalar; 182 enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime }; 183 184 Index rows = m.rows(); 185 Index cols = m.cols(); 186 187 MatrixType m1 = MatrixType::Random(rows, cols), 188 m2 = MatrixType::Random(rows, cols), 189 m3(rows, cols), 190 m4(rows, cols), 191 r1(rows, cols), 192 r2(rows, cols); 193 194 MatrixType m1up = m1.template triangularView<Upper>(); 195 MatrixType m2up = m2.template triangularView<Upper>(); 196 197 if (rows>1 && cols>1) 198 { 199 VERIFY(m1up.isUpperTriangular()); 200 VERIFY(m2up.transpose().isLowerTriangular()); 201 VERIFY(!m2.isLowerTriangular()); 202 } 203 204 // test overloaded operator+= 205 r1.setZero(); 206 r2.setZero(); 207 r1.template triangularView<Upper>() += m1; 208 r2 += m1up; 209 VERIFY_IS_APPROX(r1,r2); 210 211 // test overloaded operator= 212 m1.setZero(); 213 m1.template triangularView<Upper>() = 3 * m2; 214 m3 = 3 * m2; 215 VERIFY_IS_APPROX(m3.template triangularView<Upper>().toDenseMatrix(), m1); 216 217 218 m1.setZero(); 219 m1.template triangularView<Lower>() = 3 * m2; 220 VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1); 221 222 m1.setZero(); 223 m1.template triangularView<StrictlyUpper>() = 3 * m2; 224 VERIFY_IS_APPROX(m3.template triangularView<StrictlyUpper>().toDenseMatrix(), m1); 225 226 227 m1.setZero(); 228 m1.template triangularView<StrictlyLower>() = 3 * m2; 229 VERIFY_IS_APPROX(m3.template triangularView<StrictlyLower>().toDenseMatrix(), m1); 230 m1.setRandom(); 231 m2 = m1.template triangularView<Upper>(); 232 VERIFY(m2.isUpperTriangular()); 233 VERIFY(!m2.isLowerTriangular()); 234 m2 = m1.template triangularView<StrictlyUpper>(); 235 VERIFY(m2.isUpperTriangular()); 236 VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); 237 m2 = m1.template triangularView<UnitUpper>(); 238 VERIFY(m2.isUpperTriangular()); 239 m2.diagonal().array() -= Scalar(1); 240 VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); 241 m2 = m1.template triangularView<Lower>(); 242 VERIFY(m2.isLowerTriangular()); 243 VERIFY(!m2.isUpperTriangular()); 244 m2 = m1.template triangularView<StrictlyLower>(); 245 VERIFY(m2.isLowerTriangular()); 246 VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); 247 m2 = m1.template triangularView<UnitLower>(); 248 VERIFY(m2.isLowerTriangular()); 249 m2.diagonal().array() -= Scalar(1); 250 VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); 251 // test swap 252 m1.setOnes(); 253 m2.setZero(); 254 m2.template triangularView<Upper>().swap(m1.template triangularView<Eigen::Upper>()); 255 m3.setZero(); 256 m3.template triangularView<Upper>().setOnes(); 257 VERIFY_IS_APPROX(m2,m3); 258 } 259 260 void bug_159() 261 { 262 Matrix3d m = Matrix3d::Random().triangularView<Lower>(); 263 EIGEN_UNUSED_VARIABLE(m) 264 } 265 266 EIGEN_DECLARE_TEST(triangular) 267 { 268 int maxsize = (std::min)(EIGEN_TEST_MAX_SIZE,20); 269 for(int i = 0; i < g_repeat ; i++) 270 { 271 int r = internal::random<int>(2,maxsize); TEST_SET_BUT_UNUSED_VARIABLE(r) 272 int c = internal::random<int>(2,maxsize); TEST_SET_BUT_UNUSED_VARIABLE(c) 273 274 CALL_SUBTEST_1( triangular_square(Matrix<float, 1, 1>()) ); 275 CALL_SUBTEST_2( triangular_square(Matrix<float, 2, 2>()) ); 276 CALL_SUBTEST_3( triangular_square(Matrix3d()) ); 277 CALL_SUBTEST_4( triangular_square(Matrix<std::complex<float>,8, 8>()) ); 278 CALL_SUBTEST_5( triangular_square(MatrixXcd(r,r)) ); 279 CALL_SUBTEST_6( triangular_square(Matrix<float,Dynamic,Dynamic,RowMajor>(r, r)) ); 280 281 CALL_SUBTEST_7( triangular_rect(Matrix<float, 4, 5>()) ); 282 CALL_SUBTEST_8( triangular_rect(Matrix<double, 6, 2>()) ); 283 CALL_SUBTEST_9( triangular_rect(MatrixXcf(r, c)) ); 284 CALL_SUBTEST_5( triangular_rect(MatrixXcd(r, c)) ); 285 CALL_SUBTEST_6( triangular_rect(Matrix<float,Dynamic,Dynamic,RowMajor>(r, c)) ); 286 287 CALL_SUBTEST_100( triangular_deprecated(Matrix<float, 5, 7>()) ); 288 CALL_SUBTEST_100( triangular_deprecated(MatrixXd(r,c)) ); 289 } 290 291 CALL_SUBTEST_1( bug_159() ); 292 }