Blame ImathTest/testJacobiEigenSolver.cpp

Packit 8dc392
///////////////////////////////////////////////////////////////////////////
Packit 8dc392
//
Packit 8dc392
// Copyright (c) 2010, Industrial Light & Magic, a division of Lucas
Packit 8dc392
// Digital Ltd. LLC
Packit 8dc392
// 
Packit 8dc392
// All rights reserved.
Packit 8dc392
// 
Packit 8dc392
// Redistribution and use in source and binary forms, with or without
Packit 8dc392
// modification, are permitted provided that the following conditions are
Packit 8dc392
// met:
Packit 8dc392
// *       Redistributions of source code must retain the above copyright
Packit 8dc392
// notice, this list of conditions and the following disclaimer.
Packit 8dc392
// *       Redistributions in binary form must reproduce the above
Packit 8dc392
// copyright notice, this list of conditions and the following disclaimer
Packit 8dc392
// in the documentation and/or other materials provided with the
Packit 8dc392
// distribution.
Packit 8dc392
// *       Neither the name of Industrial Light & Magic nor the names of
Packit 8dc392
// its contributors may be used to endorse or promote products derived
Packit 8dc392
// from this software without specific prior written permission. 
Packit 8dc392
// 
Packit 8dc392
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
Packit 8dc392
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
Packit 8dc392
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
Packit 8dc392
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
Packit 8dc392
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
Packit 8dc392
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
Packit 8dc392
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
Packit 8dc392
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
Packit 8dc392
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
Packit 8dc392
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
Packit 8dc392
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
Packit 8dc392
//
Packit 8dc392
///////////////////////////////////////////////////////////////////////////
Packit 8dc392
Packit 8dc392
#include "ImathMatrix.h"
Packit 8dc392
#include "ImathMatrixAlgo.h"
Packit 8dc392
#include <iostream>
Packit 8dc392
#include <math.h>
Packit 8dc392
#include <cmath>
Packit 8dc392
#include <ctime>
Packit 8dc392
#include <cassert>
Packit 8dc392
#include <limits>
Packit 8dc392
#include <algorithm>
Packit 8dc392
Packit 8dc392
using namespace std;
Packit 8dc392
using namespace IMATH_INTERNAL_NAMESPACE;
Packit 8dc392
Packit 8dc392
const Matrix33<double> A33_1 ( 1, 0, 0, 0, 1, 0, 0, 0, 1 );
Packit 8dc392
const Matrix33<double> A33_2 ( 1, 0, 0, 0,-1, 0, 0, 0, 1 );
Packit 8dc392
const Matrix33<double> A33_3 ( 1, 0, 0, 0, 1, 0, 0, 0, 0 );
Packit 8dc392
const Matrix33<double> A33_4 ( 1, 0, 0, 0, 0, 0, 0, 0, 0 );
Packit 8dc392
const Matrix33<double> A33_5 ( 0, 0, 0, 0, 0, 0, 0, 0, 0 );
Packit 8dc392
const Matrix33<double> A33_6 ( 1, 0, 0, 0, 1e-10, 0, 0, 0, 0 );
Packit 8dc392
const Matrix33<double> A33_7 ( 1, 0, 0, 0, 1e-10, 0, 0, 0, 1e+10 );
Packit 8dc392
const Matrix33<double> A33_8 (
Packit 8dc392
     0.25058694044821,  0.49427229444416,  0.81415724537748,
Packit 8dc392
     0.49427229444416,  0.80192384710853, -0.61674948224910,
Packit 8dc392
     0.81415724537748, -0.61674948224910, -1.28486154645285);
Packit 8dc392
const Matrix33<double> A33_9 (
Packit 8dc392
     4,  -30,    60,
Packit 8dc392
   -30,  300,  -675,
Packit 8dc392
    60, -675,  1620);
Packit 8dc392
Packit 8dc392
const Matrix44<double> A44_1 ( 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1 );
Packit 8dc392
const Matrix44<double> A44_2 ( 1, 0, 0, 0, 0,-1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1 );
Packit 8dc392
const Matrix44<double> A44_3 ( 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 );
Packit 8dc392
const Matrix44<double> A44_4 ( 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 );
Packit 8dc392
const Matrix44<double> A44_5 ( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 );
Packit 8dc392
const Matrix44<double> A44_6 ( 1, 0, 0, 0, 0, 1e-20, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
Packit 8dc392
const Matrix44<double> A44_7 ( 1, 0, 0, 0, 0, 1e-20, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1e+20 );
Packit 8dc392
const Matrix44<double> A44_8 (
Packit 8dc392
     4.05747631538951,  0.16358123075600,  0.11541756047409, -1.65369223465270,
Packit 8dc392
     0.16358123075600,  0.57629829390780,  3.88542912704029,  0.92016316185369,
Packit 8dc392
     0.11541756047409,  3.88542912704029,  0.65367032943707, -0.21971103270410,
Packit 8dc392
    -1.65369223465270,  0.92016316185369, -0.21971103270410, -0.28108876552761);
Packit 8dc392
const Matrix44<double> A44_9 (
Packit 8dc392
     4,  -30,    60,   -35,
Packit 8dc392
   -30,  300,  -675,   420,
Packit 8dc392
    60, -675,  1620, -1050,
Packit 8dc392
   -35,  420, -1050,   700);
Packit 8dc392
Packit 8dc392
Packit 8dc392
template <typename TM>
Packit 8dc392
void
Packit 8dc392
verifyOrthonormal (const TM& A, const typename TM::BaseType threshold)
Packit 8dc392
{
Packit 8dc392
    const TM prod = A * A.transposed();
Packit 8dc392
    for (int i = 0; i < TM::dimensions(); ++i)
Packit 8dc392
        for (int j = 0; j < TM::dimensions(); ++j)
Packit 8dc392
            if (i == j) 
Packit 8dc392
                assert (std::abs (prod[i][j] - 1) < threshold);
Packit 8dc392
            else
Packit 8dc392
                assert (std::abs (prod[i][j]) < threshold);
Packit 8dc392
}
Packit 8dc392
Packit 8dc392
template <typename TM>
Packit 8dc392
typename TM::BaseType
Packit 8dc392
computeThreshold(const TM& A)
Packit 8dc392
{
Packit 8dc392
   typedef typename TM::BaseType T;
Packit 8dc392
   T maxAbsEntry(0);
Packit 8dc392
Packit 8dc392
   for (int i = 0; i < TM::dimensions(); ++i)
Packit 8dc392
       for (int j = 0; j < TM::dimensions(); ++j)
Packit 8dc392
           maxAbsEntry = std::max (maxAbsEntry, std::abs(A[i][j]));
Packit 8dc392
Packit 8dc392
   const T eps = std::numeric_limits<T>::epsilon();
Packit 8dc392
   maxAbsEntry = std::max(maxAbsEntry, eps);
Packit 8dc392
Packit 8dc392
   return maxAbsEntry * T(100) * eps;
Packit 8dc392
}
Packit 8dc392
Packit 8dc392
template<class TM>
Packit 8dc392
void
Packit 8dc392
testJacobiEigenSolver(const TM& A)
Packit 8dc392
{
Packit 8dc392
    using std::abs;
Packit 8dc392
Packit 8dc392
    typedef typename TM::BaseType T;
Packit 8dc392
    typedef typename TM::BaseVecType TV;
Packit 8dc392
Packit 8dc392
    const T threshold = computeThreshold(A);
Packit 8dc392
Packit 8dc392
    TM AA(A);
Packit 8dc392
    TV S;
Packit 8dc392
    TM V;
Packit 8dc392
Packit 8dc392
    jacobiEigenSolver(AA, S, V);
Packit 8dc392
Packit 8dc392
    // Orthogonality of V
Packit 8dc392
    verifyOrthonormal(V, threshold);
Packit 8dc392
 
Packit 8dc392
    // Determinant of V
Packit 8dc392
    assert(abs(V.determinant()) - 1 < threshold);
Packit 8dc392
Packit 8dc392
    // Determinant of A and S
Packit 8dc392
    TM MS;
Packit 8dc392
    for (int i = 0; i < TM::dimensions(); ++i)
Packit 8dc392
        for (int j = 0; j < TM::dimensions(); ++j)
Packit 8dc392
            if(i == j)
Packit 8dc392
                MS[i][j] = S[i];
Packit 8dc392
            else
Packit 8dc392
                MS[i][j] = 0;
Packit 8dc392
Packit 8dc392
    assert(abs(A.determinant()) - abs(MS.determinant()) <
Packit 8dc392
               threshold);
Packit 8dc392
Packit 8dc392
    // A = V * S * V^T
Packit 8dc392
    TM MA = V * MS * V.transposed();
Packit 8dc392
Packit 8dc392
    for (int i = 0; i < TM::dimensions(); ++i) 
Packit 8dc392
        for (int j =0; j < TM::dimensions(); ++j) 
Packit 8dc392
            assert(abs(A[i][j]-MA[i][j]) < threshold);
Packit 8dc392
}
Packit 8dc392
Packit 8dc392
template<class TM>
Packit 8dc392
void
Packit 8dc392
testMinMaxEigenValue(const TM& A)
Packit 8dc392
{
Packit 8dc392
  typedef typename TM::BaseVecType TV;
Packit 8dc392
  typedef typename TM::BaseType T;
Packit 8dc392
Packit 8dc392
  TV minV, maxV, S;
Packit 8dc392
  TM U, V;
Packit 8dc392
  
Packit 8dc392
  const T threshold = computeThreshold(A);
Packit 8dc392
Packit 8dc392
  {
Packit 8dc392
      TM A1(A);
Packit 8dc392
      minEigenVector(A1, minV);
Packit 8dc392
      TM A2(A);
Packit 8dc392
      maxEigenVector(A2, maxV);
Packit 8dc392
  }
Packit 8dc392
  {
Packit 8dc392
      TM A3(A);
Packit 8dc392
      jacobiSVD(A3, U, S, V);
Packit 8dc392
  }
Packit 8dc392
Packit 8dc392
  const int dim = TM::dimensions();
Packit 8dc392
Packit 8dc392
  for(int i = 0; i < dim; ++i) {
Packit 8dc392
      assert(abs(minV[i]-V[i][dim - 1]) < threshold);
Packit 8dc392
      assert(abs(maxV[i]-V[i][0]) < threshold);
Packit 8dc392
  }
Packit 8dc392
}
Packit 8dc392
Packit 8dc392
template <class T>
Packit 8dc392
void
Packit 8dc392
testJacobiTiming()
Packit 8dc392
{
Packit 8dc392
Packit 8dc392
    int rounds(100000);
Packit 8dc392
    clock_t tJacobi,tSVD, t;
Packit 8dc392
Packit 8dc392
    {
Packit 8dc392
        Matrix33<T> A,V,U;
Packit 8dc392
        Vec3<T> S;
Packit 8dc392
Packit 8dc392
        t = clock();
Packit 8dc392
        for(int i = 0; i < rounds; ++i) {
Packit 8dc392
            A = Matrix33<T>(A33_7);
Packit 8dc392
            jacobiEigenSolver(A, S, V);
Packit 8dc392
            A = Matrix33<T>(A33_8);
Packit 8dc392
            jacobiEigenSolver(A, S, V);
Packit 8dc392
        }
Packit 8dc392
        tJacobi = clock() - t;
Packit 8dc392
        cout << "Jacobi EigenSolver of 3x3 matrices took " << tJacobi << " clocks." << endl;
Packit 8dc392
Packit 8dc392
        t = clock();
Packit 8dc392
        for(int i = 0; i < rounds; ++i) {
Packit 8dc392
            A = Matrix33<T>(A33_7);
Packit 8dc392
            jacobiSVD(A, U, S, V);
Packit 8dc392
            A = Matrix33<T>(A33_8);
Packit 8dc392
            jacobiSVD(A, U, S, V);
Packit 8dc392
        }
Packit 8dc392
        tSVD = clock() - t;
Packit 8dc392
        cout << "TinySVD            of 3x3 matrices took " << tSVD << " clocks." << endl;
Packit 8dc392
        cout << (float)(tSVD-tJacobi)*100.0f/(float)(tSVD) << "% speed up." << endl;
Packit 8dc392
    }
Packit 8dc392
Packit 8dc392
    {
Packit 8dc392
        Matrix44<T> A,V,U;
Packit 8dc392
        Vec4<T> S;
Packit 8dc392
Packit 8dc392
        t = clock();
Packit 8dc392
        for(int i = 0; i < rounds; ++i) {
Packit 8dc392
            A = Matrix44<T>(A44_7);
Packit 8dc392
            jacobiEigenSolver(A, S, V);
Packit 8dc392
            A = Matrix44<T>(A44_8);
Packit 8dc392
            jacobiEigenSolver(A, S, V);
Packit 8dc392
        }
Packit 8dc392
        tJacobi = clock() - t;
Packit 8dc392
        cout << "Jacobi EigenSolver of 4x4 matrices took " << tJacobi << " clocks" << endl;
Packit 8dc392
Packit 8dc392
        t = clock();
Packit 8dc392
        for(int i = 0; i < rounds; ++i) {
Packit 8dc392
            A = Matrix44<T>(A44_7);
Packit 8dc392
            jacobiSVD(A, U, S, V);
Packit 8dc392
            A = Matrix44<T>(A44_8);
Packit 8dc392
            jacobiSVD(A, U, S, V);
Packit 8dc392
        }
Packit 8dc392
        tSVD = clock() - t;
Packit 8dc392
        cout << "TinySVD            of 4x4 matrices took " << tSVD << " clocks" << endl;
Packit 8dc392
        cout << (float)(tSVD-tJacobi)*100.0f/(float)(tSVD) << "% speed up." << endl;
Packit 8dc392
    }
Packit 8dc392
}
Packit 8dc392
Packit 8dc392
template <class T>
Packit 8dc392
void
Packit 8dc392
testJacobiEigenSolverImp()
Packit 8dc392
{
Packit 8dc392
    testJacobiEigenSolver(Matrix33<T>(A33_1));
Packit 8dc392
    testJacobiEigenSolver(Matrix33<T>(A33_2));
Packit 8dc392
    testJacobiEigenSolver(Matrix33<T>(A33_3));
Packit 8dc392
    testJacobiEigenSolver(Matrix33<T>(A33_4));
Packit 8dc392
    testJacobiEigenSolver(Matrix33<T>(A33_5));
Packit 8dc392
    testJacobiEigenSolver(Matrix33<T>(A33_6));
Packit 8dc392
    testJacobiEigenSolver(Matrix33<T>(A33_7));
Packit 8dc392
    testJacobiEigenSolver(Matrix33<T>(A33_8));
Packit 8dc392
Packit 8dc392
    testJacobiEigenSolver(Matrix44<T>(A44_1));
Packit 8dc392
    testJacobiEigenSolver(Matrix44<T>(A44_2));
Packit 8dc392
    testJacobiEigenSolver(Matrix44<T>(A44_3));
Packit 8dc392
    testJacobiEigenSolver(Matrix44<T>(A44_4));
Packit 8dc392
    testJacobiEigenSolver(Matrix44<T>(A44_5));
Packit 8dc392
    testJacobiEigenSolver(Matrix44<T>(A44_6));
Packit 8dc392
    testJacobiEigenSolver(Matrix44<T>(A44_7));
Packit 8dc392
    testJacobiEigenSolver(Matrix44<T>(A44_8));
Packit 8dc392
}
Packit 8dc392
Packit 8dc392
template <class T>
Packit 8dc392
void
Packit 8dc392
testMinMaxEigenValueImp()
Packit 8dc392
{
Packit 8dc392
    testMinMaxEigenValue(Matrix33<T>(A33_7));
Packit 8dc392
    testMinMaxEigenValue(Matrix33<T>(A33_8));
Packit 8dc392
Packit 8dc392
    testMinMaxEigenValue(Matrix44<T>(A44_7));
Packit 8dc392
    testMinMaxEigenValue(Matrix44<T>(A44_8));
Packit 8dc392
}
Packit 8dc392
Packit 8dc392
void
Packit 8dc392
testJacobiEigenSolver()
Packit 8dc392
{
Packit 8dc392
    cout << endl;
Packit 8dc392
    cout <<  "************ Testing IMATH_INTERNAL_NAMESPACE::ImathJacobiEigenSolver ************" << endl;
Packit 8dc392
    
Packit 8dc392
    cout << "Jacobi EigenSolver in single precision...";
Packit 8dc392
    testJacobiEigenSolverImp<float>();
Packit 8dc392
    cout << "PASS" << endl;
Packit 8dc392
Packit 8dc392
    cout << "Jacobi EigenSolver in double precision...";
Packit 8dc392
    testJacobiEigenSolverImp<double>();
Packit 8dc392
    cout << "PASS" << endl;
Packit 8dc392
Packit 8dc392
    cout << "Min/Max EigenValue in single precision...";
Packit 8dc392
    testMinMaxEigenValueImp<float>();
Packit 8dc392
    cout << "PASS" << endl;
Packit 8dc392
Packit 8dc392
    cout << "Min/Max EigenValue in double precision...";
Packit 8dc392
    testMinMaxEigenValueImp<double>();
Packit 8dc392
    cout << "PASS" << endl;
Packit 8dc392
Packit 8dc392
    cout << "Timing Jacobi EigenSolver in single precision...\n";
Packit 8dc392
    testJacobiTiming<float>();
Packit 8dc392
Packit 8dc392
    cout << "Timing Jacobi EigenSolver in double precision...\n";
Packit 8dc392
    testJacobiTiming<double>();
Packit 8dc392
    
Packit 8dc392
    cout << "************      ALL PASS          ************" << endl;
Packit 8dc392
}
Packit 8dc392
Packit 8dc392