1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. Eigen itself is part of the KDE project. 3 // 4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud (at) gmail.com> 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 "main.h" 11 12 template<typename MatrixType> void triangular(const MatrixType& m) 13 { 14 typedef typename MatrixType::Scalar Scalar; 15 typedef typename NumTraits<Scalar>::Real RealScalar; 16 17 RealScalar largerEps = 10*test_precision<RealScalar>(); 18 19 int rows = m.rows(); 20 int cols = m.cols(); 21 22 MatrixType m1 = MatrixType::Random(rows, cols), 23 m2 = MatrixType::Random(rows, cols), 24 m3(rows, cols), 25 m4(rows, cols), 26 r1(rows, cols), 27 r2(rows, cols); 28 29 MatrixType m1up = m1.template part<Eigen::UpperTriangular>(); 30 MatrixType m2up = m2.template part<Eigen::UpperTriangular>(); 31 32 if (rows*cols>1) 33 { 34 VERIFY(m1up.isUpperTriangular()); 35 VERIFY(m2up.transpose().isLowerTriangular()); 36 VERIFY(!m2.isLowerTriangular()); 37 } 38 39 // VERIFY_IS_APPROX(m1up.transpose() * m2, m1.upper().transpose().lower() * m2); 40 41 // test overloaded operator+= 42 r1.setZero(); 43 r2.setZero(); 44 r1.template part<Eigen::UpperTriangular>() += m1; 45 r2 += m1up; 46 VERIFY_IS_APPROX(r1,r2); 47 48 // test overloaded operator= 49 m1.setZero(); 50 m1.template part<Eigen::UpperTriangular>() = (m2.transpose() * m2).lazy(); 51 m3 = m2.transpose() * m2; 52 VERIFY_IS_APPROX(m3.template part<Eigen::LowerTriangular>().transpose(), m1); 53 54 // test overloaded operator= 55 m1.setZero(); 56 m1.template part<Eigen::LowerTriangular>() = (m2.transpose() * m2).lazy(); 57 VERIFY_IS_APPROX(m3.template part<Eigen::LowerTriangular>(), m1); 58 59 VERIFY_IS_APPROX(m3.template part<Diagonal>(), m3.diagonal().asDiagonal()); 60 61 m1 = MatrixType::Random(rows, cols); 62 for (int i=0; i<rows; ++i) 63 while (ei_abs2(m1(i,i))<1e-3) m1(i,i) = ei_random<Scalar>(); 64 65 Transpose<MatrixType> trm4(m4); 66 // test back and forward subsitution 67 m3 = m1.template part<Eigen::LowerTriangular>(); 68 VERIFY(m3.template marked<Eigen::LowerTriangular>().solveTriangular(m3).cwise().abs().isIdentity(test_precision<RealScalar>())); 69 VERIFY(m3.transpose().template marked<Eigen::UpperTriangular>() 70 .solveTriangular(m3.transpose()).cwise().abs().isIdentity(test_precision<RealScalar>())); 71 // check M * inv(L) using in place API 72 m4 = m3; 73 m3.transpose().template marked<Eigen::UpperTriangular>().solveTriangularInPlace(trm4); 74 VERIFY(m4.cwise().abs().isIdentity(test_precision<RealScalar>())); 75 76 m3 = m1.template part<Eigen::UpperTriangular>(); 77 VERIFY(m3.template marked<Eigen::UpperTriangular>().solveTriangular(m3).cwise().abs().isIdentity(test_precision<RealScalar>())); 78 VERIFY(m3.transpose().template marked<Eigen::LowerTriangular>() 79 .solveTriangular(m3.transpose()).cwise().abs().isIdentity(test_precision<RealScalar>())); 80 // check M * inv(U) using in place API 81 m4 = m3; 82 m3.transpose().template marked<Eigen::LowerTriangular>().solveTriangularInPlace(trm4); 83 VERIFY(m4.cwise().abs().isIdentity(test_precision<RealScalar>())); 84 85 m3 = m1.template part<Eigen::UpperTriangular>(); 86 VERIFY(m2.isApprox(m3 * (m3.template marked<Eigen::UpperTriangular>().solveTriangular(m2)), largerEps)); 87 m3 = m1.template part<Eigen::LowerTriangular>(); 88 VERIFY(m2.isApprox(m3 * (m3.template marked<Eigen::LowerTriangular>().solveTriangular(m2)), largerEps)); 89 90 VERIFY((m1.template part<Eigen::UpperTriangular>() * m2.template part<Eigen::UpperTriangular>()).isUpperTriangular()); 91 92 // test swap 93 m1.setOnes(); 94 m2.setZero(); 95 m2.template part<Eigen::UpperTriangular>().swap(m1); 96 m3.setZero(); 97 m3.template part<Eigen::UpperTriangular>().setOnes(); 98 VERIFY_IS_APPROX(m2,m3); 99 100 } 101 102 void selfadjoint() 103 { 104 Matrix2i m; 105 m << 1, 2, 106 3, 4; 107 108 Matrix2i m1 = Matrix2i::Zero(); 109 m1.part<SelfAdjoint>() = m; 110 Matrix2i ref1; 111 ref1 << 1, 2, 112 2, 4; 113 VERIFY(m1 == ref1); 114 115 Matrix2i m2 = Matrix2i::Zero(); 116 m2.part<SelfAdjoint>() = m.part<UpperTriangular>(); 117 Matrix2i ref2; 118 ref2 << 1, 2, 119 2, 4; 120 VERIFY(m2 == ref2); 121 122 Matrix2i m3 = Matrix2i::Zero(); 123 m3.part<SelfAdjoint>() = m.part<LowerTriangular>(); 124 Matrix2i ref3; 125 ref3 << 1, 0, 126 0, 4; 127 VERIFY(m3 == ref3); 128 129 // example inspired from bug 159 130 int array[] = {1, 2, 3, 4}; 131 Matrix2i::Map(array).part<SelfAdjoint>() = Matrix2i::Random().part<LowerTriangular>(); 132 133 std::cout << "hello\n" << array << std::endl; 134 } 135 136 void test_eigen2_triangular() 137 { 138 CALL_SUBTEST_8( selfadjoint() ); 139 for(int i = 0; i < g_repeat ; i++) { 140 CALL_SUBTEST_1( triangular(Matrix<float, 1, 1>()) ); 141 CALL_SUBTEST_2( triangular(Matrix<float, 2, 2>()) ); 142 CALL_SUBTEST_3( triangular(Matrix3d()) ); 143 CALL_SUBTEST_4( triangular(MatrixXcf(4, 4)) ); 144 CALL_SUBTEST_5( triangular(Matrix<std::complex<float>,8, 8>()) ); 145 CALL_SUBTEST_6( triangular(MatrixXd(17,17)) ); 146 CALL_SUBTEST_7( triangular(Matrix<float,Dynamic,Dynamic,RowMajor>(5, 5)) ); 147 } 148 } 149
