[vtkusers] Jacobi eigenvectors incorrect?

[Bryn Lloyd]

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;
}