[vtkusers] Jacobi eigenvectors incorrect?

[Bill Lorensen]

According to the documentation in vtkMath.cxx, the code
vtkJacobiN finds the solution of an nxn real SYMMETRIC matrix.
Can you try your tests with a symmetric matrix and compare the results
between vtk and matlab.

Bill

On Tue, Oct 27, 2009 at 7:56 AM, David Doria <daviddoria+vtk at gmail.com> wrote:
> 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
> _______________________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the VTK FAQ at: http://www.vtk.org/Wiki/VTK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.vtk.org/mailman/listinfo/vtkusers
>