From e0ea25fc213a49198190ead2cdc9da3d3b59f21e Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 1 Sep 2010 12:59:38 +0200 Subject: add missing copyrights --- bench/eig33.cpp | 168 ++++++++++++++++++++++++++++++++++++-------------------- 1 file changed, 107 insertions(+), 61 deletions(-) (limited to 'bench/eig33.cpp') diff --git a/bench/eig33.cpp b/bench/eig33.cpp index 2016c2c01..df07ad79d 100644 --- a/bench/eig33.cpp +++ b/bench/eig33.cpp @@ -1,3 +1,56 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +// The computeRoots function included in this is based on materials +// covered by the following copyright and license: +// +// Geometric Tools, LLC +// Copyright (c) 1998-2010 +// Distributed under the Boost Software License, Version 1.0. +// +// Permission is hereby granted, free of charge, to any person or organization +// obtaining a copy of the software and accompanying documentation covered by +// this license (the "Software") to use, reproduce, display, distribute, +// execute, and transmit the Software, and to prepare derivative works of the +// Software, and to permit third-parties to whom the Software is furnished to +// do so, all subject to the following: +// +// The copyright notices in the Software and this entire statement, including +// the above license grant, this restriction and the following disclaimer, +// must be included in all copies of the Software, in whole or in part, and +// all derivative works of the Software, unless such copies or derivative +// works are solely in the form of machine-executable object code generated by +// a source language processor. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT +// SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE +// FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, +// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. + #include #include #include @@ -8,56 +61,49 @@ using namespace Eigen; using namespace std; template -inline void computeRoots (const Matrix& rkA, Roots& adRoot) +inline void computeRoots(const Matrix& m, Roots& roots) { typedef typename Matrix::Scalar Scalar; - const Scalar msInv3 = 1.0/3.0; - const Scalar msRoot3 = ei_sqrt(Scalar(3.0)); - - Scalar dA00 = rkA(0,0); - Scalar dA01 = rkA(0,1); - Scalar dA02 = rkA(0,2); - Scalar dA11 = rkA(1,1); - Scalar dA12 = rkA(1,2); - Scalar dA22 = rkA(2,2); + const Scalar s_inv3 = 1.0/3.0; + const Scalar s_sqrt3 = ei_sqrt(Scalar(3.0)); // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0. The // eigenvalues are the roots to this equation, all guaranteed to be // real-valued, because the matrix is symmetric. - Scalar dC0 = dA00*dA11*dA22 + Scalar(2)*dA01*dA02*dA12 - dA00*dA12*dA12 - dA11*dA02*dA02 - dA22*dA01*dA01; - Scalar dC1 = dA00*dA11 - dA01*dA01 + dA00*dA22 - dA02*dA02 + dA11*dA22 - dA12*dA12; - Scalar dC2 = dA00 + dA11 + dA22; + 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); + 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); + Scalar c2 = m(0,0) + m(1,1) + m(2,2); // Construct the parameters used in classifying the roots of the equation // and in solving the equation for the roots in closed form. - Scalar dC2Div3 = dC2*msInv3; - Scalar dADiv3 = (dC1 - dC2*dC2Div3)*msInv3; - if (dADiv3 > Scalar(0)) - dADiv3 = Scalar(0); + Scalar c2_over_3 = c2*s_inv3; + Scalar a_over_3 = (c1 - c2*c2_over_3)*s_inv3; + if (a_over_3 > Scalar(0)) + a_over_3 = Scalar(0); - Scalar dMBDiv2 = Scalar(0.5)*(dC0 + dC2Div3*(Scalar(2)*dC2Div3*dC2Div3 - dC1)); + Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1)); - Scalar dQ = dMBDiv2*dMBDiv2 + dADiv3*dADiv3*dADiv3; - if (dQ > Scalar(0)) - dQ = Scalar(0); + Scalar q = half_b*half_b + a_over_3*a_over_3*a_over_3; + if (q > Scalar(0)) + q = Scalar(0); // Compute the eigenvalues by solving for the roots of the polynomial. - Scalar dMagnitude = ei_sqrt(-dADiv3); - Scalar dAngle = std::atan2(ei_sqrt(-dQ),dMBDiv2)*msInv3; - Scalar dCos = ei_cos(dAngle); - Scalar dSin = ei_sin(dAngle); - adRoot(0) = dC2Div3 + 2.f*dMagnitude*dCos; - adRoot(1) = dC2Div3 - dMagnitude*(dCos + msRoot3*dSin); - adRoot(2) = dC2Div3 - dMagnitude*(dCos - msRoot3*dSin); + Scalar rho = ei_sqrt(-a_over_3); + Scalar theta = std::atan2(ei_sqrt(-q),half_b)*s_inv3; + Scalar cos_theta = ei_cos(theta); + Scalar sin_theta = ei_sin(theta); + roots(0) = c2_over_3 + Scalar(2)*rho*cos_theta; + roots(1) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta); + roots(2) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta); // Sort in increasing order. - if (adRoot(0) >= adRoot(1)) - std::swap(adRoot(0),adRoot(1)); - if (adRoot(1) >= adRoot(2)) + if (roots(0) >= roots(1)) + std::swap(roots(0),roots(1)); + if (roots(1) >= roots(2)) { - std::swap(adRoot(1),adRoot(2)); - if (adRoot(0) >= adRoot(1)) - std::swap(adRoot(0),adRoot(1)); + std::swap(roots(1),roots(2)); + if (roots(0) >= roots(1)) + std::swap(roots(0),roots(1)); } } @@ -65,21 +111,21 @@ template void eigen33(const Matrix& mat, Matrix& evecs, Vector& evals) { typedef typename Matrix::Scalar Scalar; - // Scale the matrix so its entries are in [-1,1]. The scaling is applied - // only when at least one matrix entry has magnitude larger than 1. + // Scale the matrix so its entries are in [-1,1]. The scaling is applied + // only when at least one matrix entry has magnitude larger than 1. - Scalar scale = mat.cwiseAbs()/*.template triangularView()*/.maxCoeff(); - scale = std::max(scale,Scalar(1)); - Matrix scaledMat = mat / scale; + Scalar scale = mat.cwiseAbs()/*.template triangularView()*/.maxCoeff(); + scale = std::max(scale,Scalar(1)); + Matrix scaledMat = mat / scale; - // Compute the eigenvalues -// scaledMat.setZero(); - computeRoots(scaledMat,evals); + // Compute the eigenvalues +// scaledMat.setZero(); + computeRoots(scaledMat,evals); - // compute the eigen vectors - // here we assume 3 differents eigenvalues + // compute the eigen vectors + // **here we assume 3 differents eigenvalues** - // "optimized version" which appears to be slower with gcc! + // "optimized version" which appears to be slower with gcc! // Vector base; // Scalar alpha, beta; // base << scaledMat(1,0) * scaledMat(2,1), @@ -93,22 +139,22 @@ void eigen33(const Matrix& mat, Matrix& evecs, Vector& evals) // } // evecs.col(2) = evecs.col(0).cross(evecs.col(1)).normalized(); - // naive version - Matrix tmp; - tmp = scaledMat; - tmp.diagonal().array() -= evals(0); - evecs.col(0) = tmp.row(0).cross(tmp.row(1)).normalized(); - - tmp = scaledMat; - tmp.diagonal().array() -= evals(1); - evecs.col(1) = tmp.row(0).cross(tmp.row(1)).normalized(); - - tmp = scaledMat; - tmp.diagonal().array() -= evals(2); - evecs.col(2) = tmp.row(0).cross(tmp.row(1)).normalized(); - - // Rescale back to the original size. - evals *= scale; + // naive version + Matrix tmp; + tmp = scaledMat; + tmp.diagonal().array() -= evals(0); + evecs.col(0) = tmp.row(0).cross(tmp.row(1)).normalized(); + + tmp = scaledMat; + tmp.diagonal().array() -= evals(1); + evecs.col(1) = tmp.row(0).cross(tmp.row(1)).normalized(); + + tmp = scaledMat; + tmp.diagonal().array() -= evals(2); + evecs.col(2) = tmp.row(0).cross(tmp.row(1)).normalized(); + + // Rescale back to the original size. + evals *= scale; } int main() -- cgit v1.2.3