Re: [ublas] [bindings][traits] is_vector, is_matrix

Subject: Re: [ublas] [bindings][traits] is_vector, is_matrix
From: Rutger ter Borg (rutger_at_[hidden])
Date: 2009-03-12 03:13:46

Jesse Perla wrote:

> Speaking of which...  The bindings for ?tritri don't seem to work as one
> would expect for triangular matrices.  The following code works:
> ublas::matrix<double, ublas::column_major> (A2,2);
> A(0,0) = 1.0;
> A(1,0) = 1.0; A(1, 1) = 2.0;
> A(0, 1) = 0.0; //Shouldn't need, but good to set
> lapack::trtri('L', 'N', A);

This is right, trtri is meant for normal dense matrices (=storage format) of
which you say they are triangular. I've added support for uplo detection in
trtri, to use it now please use ublas' triangular adaptor to indicate
whether it is upper/lower. 

I'm not sure whether it is possible to do the diag detection through the
traits at the moment -- Karl?
 
> But the following won't compile, which has the triangular matrix one would
> expect to work:
> 
> ublas::triangular_matrix<double, ublas::lower, ublas::column_major>
> A(2,2); A(0,0) = 1.0;
> A(1,0) = 1.0; A(1, 1) = 2.0;
> lapack::trtri('L', 'N', A);

I think ublas uses a packed storage format for triangular matrices. Please
try tptri, I've just added it to the computational section. However, I
think the traits do not support triangular matrices yet?

> The new bindings seem to compile and my basic test cases worked.  I still
> don't really know what is going on with the workspace stuff, but should
> the following code (which compiles) work? (it pops out a different number
> than matlab, but I could be misinterpreting the results):
> 
> ublas::matrix<double, ublas::column_major> A(2,2);
> //First need to call getrf()
> std::vector<int> ipiv (2);  // pivot vector
> lapack::getrf(A, ipiv);
> //Then can get the inverse
> lapack::getri(A, ipiv, lapack::optimal_workspace() );
> 
> ... optimal_workspace() is compiling now, but is it working?

If you would like to see what workspace size has been chosen, you could
insert a 

std::cout << "Workspace real size: " << traits::detail::to_int(
opt_size_work ) << std::endl;

at line 110/111 of getri.

Cheers,

Rutger