1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr> 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// The computeRoots function included in this is based on materials 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// covered by the following copyright and license: 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Geometric Tools, LLC 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (c) 1998-2010 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Distributed under the Boost Software License, Version 1.0. 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Permission is hereby granted, free of charge, to any person or organization 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// obtaining a copy of the software and accompanying documentation covered by 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// this license (the "Software") to use, reproduce, display, distribute, 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// execute, and transmit the Software, and to prepare derivative works of the 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Software, and to permit third-parties to whom the Software is furnished to 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// do so, all subject to the following: 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// The copyright notices in the Software and this entire statement, including 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// the above license grant, this restriction and the following disclaimer, 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// must be included in all copies of the Software, in whole or in part, and 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// all derivative works of the Software, unless such copies or derivative 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// works are solely in the form of machine-executable object code generated by 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// a source language processor. 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// DEALINGS IN THE SOFTWARE. 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <iostream> 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <Eigen/Core> 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <Eigen/Eigenvalues> 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <Eigen/Geometry> 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <bench/BenchTimer.h> 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathusing namespace Eigen; 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathusing namespace std; 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Matrix, typename Roots> 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline void computeRoots(const Matrix& m, Roots& roots) 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Matrix::Scalar Scalar; 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar s_inv3 = 1.0/3.0; 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar s_sqrt3 = internal::sqrt(Scalar(3.0)); 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0. The 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // eigenvalues are the roots to this equation, all guaranteed to be 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // real-valued, because the matrix is symmetric. 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar c0 = m(0,0)*m(1,1)*m(2,2) + Scalar(2)*m(0,1)*m(0,2)*m(1,2) - m(0,0)*m(1,2)*m(1,2) - m(1,1)*m(0,2)*m(0,2) - m(2,2)*m(0,1)*m(0,1); 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar c1 = m(0,0)*m(1,1) - m(0,1)*m(0,1) + m(0,0)*m(2,2) - m(0,2)*m(0,2) + m(1,1)*m(2,2) - m(1,2)*m(1,2); 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar c2 = m(0,0) + m(1,1) + m(2,2); 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Construct the parameters used in classifying the roots of the equation 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // and in solving the equation for the roots in closed form. 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar c2_over_3 = c2*s_inv3; 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar a_over_3 = (c1 - c2*c2_over_3)*s_inv3; 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (a_over_3 > Scalar(0)) 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath a_over_3 = Scalar(0); 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1)); 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar q = half_b*half_b + a_over_3*a_over_3*a_over_3; 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (q > Scalar(0)) 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath q = Scalar(0); 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute the eigenvalues by solving for the roots of the polynomial. 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar rho = internal::sqrt(-a_over_3); 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar theta = std::atan2(internal::sqrt(-q),half_b)*s_inv3; 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar cos_theta = internal::cos(theta); 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar sin_theta = internal::sin(theta); 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath roots(0) = c2_over_3 + Scalar(2)*rho*cos_theta; 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath roots(1) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta); 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath roots(2) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta); 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Sort in increasing order. 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (roots(0) >= roots(1)) 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::swap(roots(0),roots(1)); 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (roots(1) >= roots(2)) 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::swap(roots(1),roots(2)); 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (roots(0) >= roots(1)) 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::swap(roots(0),roots(1)); 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Matrix, typename Vector> 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid eigen33(const Matrix& mat, Matrix& evecs, Vector& evals) 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Matrix::Scalar Scalar; 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Scale the matrix so its entries are in [-1,1]. The scaling is applied 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // only when at least one matrix entry has magnitude larger than 1. 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar scale = mat.cwiseAbs()/*.template triangularView<Lower>()*/.maxCoeff(); 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath scale = std::max(scale,Scalar(1)); 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix scaledMat = mat / scale; 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute the eigenvalues 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// scaledMat.setZero(); 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath computeRoots(scaledMat,evals); 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // compute the eigen vectors 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // **here we assume 3 differents eigenvalues** 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // "optimized version" which appears to be slower with gcc! 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Vector base; 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Scalar alpha, beta; 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// base << scaledMat(1,0) * scaledMat(2,1), 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// scaledMat(1,0) * scaledMat(2,0), 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// -scaledMat(1,0) * scaledMat(1,0); 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for(int k=0; k<2; ++k) 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// { 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// alpha = scaledMat(0,0) - evals(k); 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// beta = scaledMat(1,1) - evals(k); 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// evecs.col(k) = (base + Vector(-beta*scaledMat(2,0), -alpha*scaledMat(2,1), alpha*beta)).normalized(); 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// } 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// evecs.col(2) = evecs.col(0).cross(evecs.col(1)).normalized(); 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// // naive version 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Matrix tmp; 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// tmp = scaledMat; 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// tmp.diagonal().array() -= evals(0); 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// evecs.col(0) = tmp.row(0).cross(tmp.row(1)).normalized(); 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// tmp = scaledMat; 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// tmp.diagonal().array() -= evals(1); 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// evecs.col(1) = tmp.row(0).cross(tmp.row(1)).normalized(); 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// tmp = scaledMat; 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// tmp.diagonal().array() -= evals(2); 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// evecs.col(2) = tmp.row(0).cross(tmp.row(1)).normalized(); 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // a more stable version: 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if((evals(2)-evals(0))<=Eigen::NumTraits<Scalar>::epsilon()) 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath evecs.setIdentity(); 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix tmp; 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp = scaledMat; 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp.diagonal ().array () -= evals (2); 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath evecs.col (2) = tmp.row (0).cross (tmp.row (1)).normalized (); 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp = scaledMat; 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp.diagonal ().array () -= evals (1); 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath evecs.col(1) = tmp.row (0).cross(tmp.row (1)); 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar n1 = evecs.col(1).norm(); 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(n1<=Eigen::NumTraits<Scalar>::epsilon()) 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath evecs.col(1) = evecs.col(2).unitOrthogonal(); 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath evecs.col(1) /= n1; 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // make sure that evecs[1] is orthogonal to evecs[2] 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath evecs.col(1) = evecs.col(2).cross(evecs.col(1).cross(evecs.col(2))).normalized(); 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath evecs.col(0) = evecs.col(2).cross(evecs.col(1)); 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Rescale back to the original size. 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath evals *= scale; 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathint main() 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BenchTimer t; 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int tries = 10; 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int rep = 400000; 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix3f Mat; 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Vector3f Vec; 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Mat A = Mat::Random(3,3); 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath A = A.adjoint() * A; 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SelfAdjointEigenSolver<Mat> eig(A); 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BENCH(t, tries, rep, eig.compute(A)); 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << "Eigen: " << t.best() << "s\n"; 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Mat evecs; 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Vec evals; 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BENCH(t, tries, rep, eigen33(A,evecs,evals)); 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << "Direct: " << t.best() << "s\n\n"; 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cerr << "Eigenvalue/eigenvector diffs:\n"; 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cerr << (evals - eig.eigenvalues()).transpose() << "\n"; 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int k=0;k<3;++k) 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(evecs.col(k).dot(eig.eigenvectors().col(k))<0) 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath evecs.col(k) = -evecs.col(k); 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cerr << evecs - eig.eigenvectors() << "\n\n"; 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 197