Re: LU decomposition of complex dense matrices

Thanks Gunter, for setting me straight about indices
starting from zero and about max_index.  

Below is a sample program which demonstrates
lu_factorize working on a 2x2 complex matrix.  I
believe that the output is incorrect.  How can I fix
it?

Dave



// -*- c++ -*-
//
// $COPYRIGHT$
//
//===========================================================================
#include "C:/Include/mtl/workaround.h"
#include "C:/Include/mtl/utils.h" 
#include "C:/Include/mtl/linalg_vec.h"
#include "C:/Include/mtl/mtl.h"
#include "C:/Include/mtl/lu.h"
using namespace mtl;
using namespace std;
typedef external_vec<double> Vec;

int main()
{
	typedef matrix<complex<double>, rectangle<>,dense<>,
row_major>::type MyMatrix;
	const int N = 2;
	MyMatrix A(N, N);
	A(0,0) = complex<double>(1);
	A(0,1) = complex<double>(0,1);
	A(1,0) = complex<double>(0,-1);
	A(1,1) = complex<double>(-1);
	int zero_vec[] = {0,0};
	external_vec<int> ipvt(zero_vec,N);
/*  Sample Output
		Input Matrix to be factorized:
		2x2
		[
		[(1,0),(0,1)],
		[(0,-1),(-1,0)]
		]
		Factorization:
		2x2
		[
		[(1,0),(0,1)],
		[(0,-1),(0,0)]
		]
		Press any key to continue
*/  
/* Desired Output 
		Input Matrix to be factorized:
		2x2
		[
		[(1,0),(0,1)],
		[(0,-1),(-1,0)]
		]
		Factorization:
		2x2
		[
		[(1,0),(0,1)],
		[(0,-1),(-2,0)]
		]
		Press any key to continue
*/  
	cout << "Input Matrix to be factorized:  " << endl;
	print_all_matrix(A);
	lu_factorize(A, ipvt);
	cout << "Factorization:  " << endl;
	print_all_matrix(A);
  return 0;
}
 

=====
[email protected]
David Lidsky


__________________________________________________
Do you Yahoo!?
New DSL Internet Access from SBC & Yahoo!
http://sbc.yahoo.com
_______________________________________________
This list is archived at http://www.osl.iu.edu/MailArchives/mtl-devel/