[Mplapack-devel] possibly a bug in Cheev

[Akira S. <aki...@ni...>]

Dear Dr. Maho NAKATA,

I am a postdoc working on quantum information and the matrices
I have to handle are quite often Hermitian matrices with some
eigenvalues degenerate or close to each other. I am trying several
libraries as well as writing my own library for matrix calculations
required in my research.

I would like to report some issue as I think this is possibly a
bug of MPACK. Or, it might be a problem in my code or the environment.

I tried Cheev of MPACK with GMP but could not get correct results
for some simple matrices. For example, it does not diagonalize the
3 x 3 matrix

1 0 1
0 0 0
1 0 1

and gives me "#Rlaev2 Checkpoint 13 Not checked". It works fine
for diag(1,0,0) and a matrix filled with 1's among others.

I am using mpack-0.7.0 on Fedora 15 where GMP 4.3.2 exists.

Thanks for your attentions.
Regards,
Akira SAITOH, NII, Japan

===========================================
[saitoh@localhost test_cheev]$ g++ test.cpp -L/usr/local/lib 
-lmlapack_gmp -lmblas_gmp -lmblas_gmp_ref -lmlapack_gmp_ref -lgmp -lmpfr 
-lgmpxx -lgomp -g
/usr/bin/ld: warning: libgmp.so.3, needed by 
/usr/lib/gcc/x86_64-redhat-linux/4.6.3/../../../../lib64/libmpfr.so, may 
conflict with libgmp.so.10
[saitoh@localhost test_cheev]$ LD_LIBRARY_PATH=/usr/local/lib ./a.out 1 
0 1 0 0 0 1 0 1
A =
1 0 1
0 0 0
1 0 1
#Rlaev2 Checkpoint 13 Not checked
[saitoh@localhost test_cheev]$ LD_LIBRARY_PATH=/usr/local/lib ./a.out 1 
0 0 0 0 0 0 0 0
A =
1 0 0
0 0 0
0 0 0
Eigenvalues:
0 0 1
Unitary matrix:
0+i*0 0+i*0 1+i*0
1+i*0 0+i*0 0+i*0
0+i*0 1+i*0 0+i*0
===========================================

test.cpp:

===========================================
#include <gmpxx.h>
#include <mpack/mpack_config.h>
#include <mpack/mpc_class.h>
#include <mpack/mlapack_gmp.h>
#include <cstdlib>
#include <iostream>
int main (int argc, char* argv[])
{
   mpf_set_default_prec(512);
   mpackint n = 3;
   mpc_class *A = new mpc_class [n * n];
   mpackint lda = n;
   mpf_class *w = new mpf_class [n];
   mpc_class *work = new mpc_class [2 * n - 1];
   mpackint lwork = 2 * n - 1;
   mpf_class *rwork = new mpf_class [3 * n - 2];
   mpackint info;

   for (int i = 1; i < argc && i < n * n + 1; i++)
     A[i-1] = ::atof(argv[i]);

   std::cout << "A = " << std::endl;
   for (int i = 0; i < n; i ++)
   {
     for (int j = 0; j < n; j++)
       std::cout << A[i + n * j].real().get_d() << " ";
     std::cout << std::endl;
   }

   Cheev ("V", "U", n, A, lda, w, work, lwork, rwork, &info);

   std::cout << "Eigenvalues: " << std::endl;
   for (int i = 0; i < n; i++)
     std::cout << w[i].get_d() << " ";
   std::cout << std::endl;

   std::cout << "Unitary matrix:" << std::endl;
   for (int i = 0; i < n; i ++)
   {
     for (int j = 0; j < n; j++)
       std::cout << A[i + n * j].real().get_d() << "+i*"
                 << A[i + n * j].imag().get_d() << " ";
     std::cout << std::endl;
   }

   delete [] A, w, work, rwork;
   return 0;
}
===========================================