lu for complex<double> bug?

Mateusz Koreń <[email protected]> Mon, 08 May 2006 11:38:55 +0200
Newsgroups gmane.comp.lib.mtl.devel
Organization Nowa Poczta Wirtualnej Polski S.A. http://www.wp.pl/
Message-ID <[email protected]>
I'm trying to use mtl as a inverter for my tvmet matrices and i found that something gets wrong when matrix elements have imaginary parts. 
I checked how it works with lu_factorization.cc from your website and it seems to give something like this (sorry for an inelegant code):

#include <iostream>
#include <complex>
using namespace std;

#include <mtl/lu.h>
#include <mtl/matrix.h>
#include <mtl/dense1D.h>
#include <mtl/utils.h>
#include <mtl/blais.h>

int main() {
  using namespace mtl;

  typedef complex<double> comp;
  typedef matrix<comp,
                 rectangle<>,
                 dense<external>,
                 column_major>::type Matrix;
  const int N = 3;
  comp da [] = {3,comp(1,1),0,0,3,0,comp(2,0),0,3};
  comp da2 [] = {3,comp(1,1),0,0,3,0,comp(2,0),0,3};
  comp db [] = {3,comp(1,1),0,0,3,0,comp(2,2),0,3};
  comp db2 [] = {3,comp(1,1),0,0,3,0,comp(2,2),0,3};
  comp resulta[9];
  comp resultb[9];
  comp onea[9];
  comp oneb[9];
  Matrix A2(da2,N,N);
  Matrix A(da,N,N);
  Matrix B2(db2,N,N);
  Matrix B(db,N,N);
  Matrix A_inv(resulta,N,N);
  Matrix B_inv(resultb,N,N);
  Matrix identa(onea,N,N);
  Matrix identb(oneb,N,N);
  dense1D<int> pivotsa(N,0);
  dense1D<int> pivotsb(N,0);

  lu_factor(A, pivotsa);
  lu_inverse(A,pivotsa,A_inv);

  lu_factor(B, pivotsb);
  lu_inverse(B,pivotsb,B_inv);

  mult(A2,A_inv,identa);
  mult(B2,B_inv,identb);
  cout<<"A*A^-1 (should be identity):"<<endl;
  print_all_matrix(identa);
  cout<<"B*B^-1 (should be identity):"<<endl;
  print_all_matrix(identb);
  return 0;
}

resulting in:

 A*A^-1 (should be identity):
3x3
[
[(1,0),(0,0),(0,0)],
[(0,0),(1,0),(0,0)],
[(0,0),(0,0),(1,0)]
]
B*B^-1 (should be identity):
3x3
[
[(1,0),(0,0),(0,0)],
[(0,0),(1,0),(0.444444,-0.444444)],
[(0,0),(0,0),(1,0)]
]

and so on with more "imaginary" matrices. Seems like a bug to me. Or maybe I'm doing something wrong? I've got Debian/GNU Linux 3.1 with gcc 4.0

Best regards,
Mateusz

PS I also checked the internal storage way - with the same result.

----------------------------------------------------
Koncert zespołu TOOL! 
24 czerwca w katowickim Spodku!
http://klik.wp.pl/?adr=http%3A%2F%2Fadv.reklama.wp.pl%2Fas%2Ftool.html&sid=742


_______________________________________________
This list is archived at http://www.osl.iu.edu/MailArchives/mtl-devel/