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