text: // This is core/vnl/algo/tests/test_generalized_eigensystem.cxx
  #include <iostream>
  #include "testlib/testlib_test.h"
  //:
  // \file
  // \brief test program for generalized eigensystem routines.
  // \author Andrew W. Fitzgibbon, Oxford RRG.
  // \date 29 Aug 96
  
  //-----------------------------------------------------------------------------
  
  #include <vnl/algo/vnl_generalized_eigensystem.h>
  
  void test_generalized_eigensystem()
  {
    double Sdata[36] = {
     30.0000,   -3.4273,   13.9254,   13.7049,   -2.4446,   20.2380,
     -3.4273,   13.7049,   -2.4446,    1.3659,    3.6702,   -0.2282,
     13.9254,   -2.4446,   20.2380,    3.6702,   -0.2282,   28.6779,
     13.7049,    1.3659,    3.6702,   12.5273,   -1.6045,    3.9419,
     -2.4446,    3.6702,   -0.2282,   -1.6045,    3.9419,    2.5821,
     20.2380,   -0.2282,   28.6779,    3.9419,    2.5821,   44.0636,
    };
    double Cdata[36] = {
      0,  0,  0,  0,  0,  0,
      0,  0,  0,  0,  0,  0,
      0,  0,  0,  0,  0,  0,
      0,  0,  0,  0,  0,  2,
      0,  0,  0,  0, -1,  0,
      0,  0,  0,  2,  0,  0,
    };
  
    vnl_matrix<double> S(Sdata, 6,6);
    vnl_matrix<double> C(Cdata, 6,6);
  
    vnl_generalized_eigensystem gev(C, S);
  
    std::cout << "V = " << gev.V << std::endl
             << "D = " << gev.D << std::endl
             << "residual = " << C * gev.V - S * gev.V * gev.D << std::endl;
    double err = (C * gev.V - S * gev.V * gev.D).fro_norm();
    std::cout << "Recomposition residual = " << err << std::endl;
  
    TEST_NEAR("Recomposition residual < 1e-12", err, 0.0, 1e-12);
  }
  
  TESTMAIN(test_generalized_eigensystem);
  