[vtkusers] Jacobi eigenvectors incorrect?

[David Doria]

On Tue, Oct 27, 2009 at 4:18 AM, Bryn Lloyd <blloyd at vision.ee.ethz.ch> wrote:
> Hi David,
>
> I agree. The result is different to what I get in Matlab for a different
> matrix. vtkMath::JacobiN returns the same result as vtkMath::Jacobi.
>
>
> My test programm:
>
>
>
>
> #include <vtkMath.h>
>
>
> template<class TReal>
> TReal **create_matrix(long nrow, long ncol) {
>  typedef TReal* TRealPointer;
>  TReal **m = new TRealPointer[nrow];
>
>  TReal* block = (TReal*)calloc(nrow*ncol, sizeof(TReal));
>  m[0] = block;
>  for ( int row = 1; row < nrow; ++row ) {
>    m[ row ] = &block[ row * ncol ];
>  }
>  return m;
> }
> /* free a TReal matrix allocated with matrix() */
> template<class TReal>
> void free_matrix(TReal **m) {
>  free (m[0]);
>  delete[] m;
> }
>
>
> int main(int argc, char ** argv)
> {
>  // Matrix
>  // 1  3  2
>  // 1 -2  3
>  // 1  2 -3
>
>  double **a = create_matrix<double>(3,3);
>  a[0][0] = 1; a[0][1] = 3;  a[0][2] = 2;
>  a[1][0] = 1; a[1][1] = -2; a[1][2] = 3;
>  a[2][0] = 1; a[2][1] = 2;  a[2][2] = -3;
>
>  double **w = create_matrix<double>(3,3);
>  double v[3];
>
>  printf("Jacobi\n");
>  vtkMath::Jacobi(a,v,w);
>  for(int i=0; i<3; i++)
>    printf("%g \t%g %g %g\n",v[i],w[i][0],w[i][1],w[i][2]);
>
>
>  a[0][0] = 1; a[0][1] = 3;  a[0][2] = 2;
>  a[1][0] = 1; a[1][1] = -2; a[1][2] = 3;
>  a[2][0] = 1; a[2][1] = 2;  a[2][2] = -3;
>
>  printf("\nJacobiN\n");
>  vtkMath::JacobiN(a,3,v,w);
>  for(int i=0; i<3; i++)
>    printf("%g \t%g %g %g\n",v[i],w[i][0],w[i][1],w[i][2]);
>
> /*
> Output from this program:
> Jacobi
> 4.369           0.729826 -0.677467 0.0916155
> -2.78956        0.540772 0.490118 -0.68363
> -5.57944        0.418234 0.548474 0.724056
>
> JacobiN
> 4.369           0.729826 -0.677467 0.0916155
> -2.78956        0.540772 0.490118 -0.68363
> -5.57944        0.418234 0.548474 0.724056
> */
>
>
> /*
> % Matlab result:
>
> a = [[1,3,2];[1,-2,3];[1,2,-3]];
> [v,w] = eig(a)
>
> v =
>
>    0.8924    0.7816    0.1413
>    0.3570   -0.5541   -0.7279
>    0.2761   -0.2864    0.6709
>
>
> w =
>
>    2.8190         0         0
>         0   -1.8598         0
>         0         0   -4.9593
> */
>  return 0;
> }

Thanks for taking a look Bryn. JacobiN seems like it is just a
particular iteration of Jacobi()? Just curious why you would choose to
compare Jacobi to JacobiN?

This seems like a mildly big problem... as any algorithms that are
using eigenvalues/vectors could be producing bogus results silently.
Is anyone familiar with the internals of the Jacobi functions?

Thanks,

David