|
|
Ublas : |
Subject: [ublas] [numeric_bindings] Call LAPACK and BLAS with C arrays
From: sguazt (marco.guazzone_at_[hidden])
Date: 2012-01-30 16:36:02
Hello,
Is it possible to call any of the ::boost::numeric::bindings::lapack
or ::boost::numeric::bindings::blas functions by using plain C arrays?
In my case a C array can represent either a vector or a linearized
matrix (in column-major order).
For instance, take the LAPACK ?POSV subroutine.
Currently I use this "trick":
namespace ublas = ::boost::numeric::ublas;
namespace bindings = ::boost::numeric::bindings;
typedef double real_type;
typedef ublas::array_adaptor<real_type> array_type;
typedef ublas::matrix<real_type,ublas::column_major,array_type> matrix_type;
typedef ublas::vector<real_type,array_type> vector_type;
matrix_type AA(n, n, array_type(n*n, A));
matrix_type AA(n, m, array_type(n*m, B));
// Call DPOSV("L", &n, &n, A, &n, B, &m, &info);
bindings::lapack::posv(bindings::lower(AA), BB);
This trick works in part. The problem is in A which seems not to be updated.
That is, if I print AA the result is OK, but if I print A the result is not OK.
Here below is the output I get by running the sample code attached to
this email.
Note the difference between A and AA in the "After POSV" prints.
The correct result is that of AA.
Before POSV -- A=[6,6]((2.0013,0.0000,0.0000,0.0000,0.0000,0.0000),(0.0000,0.0000,0.0000,0.0000,2.0013,0.0000),(0.0000,0.0000,0.0000,2.0013,0.0000,0.0000),(0.0000,0.0000,2.0013,0.0000,0.0000,0.0000),(0.0000,2.0013,0.0000,0.0000,0.0000,0.0000),(2.0013,0.0000,0.0000,0.0000,0.0000,0.0000))
Before POSV -- B=[6,3]((0.7627,0.1149,0.0025),(0.1149,0.7652,0.1149),(0.0025,0.1149,0.7652),(0.0000,0.0025,0.1149),(0.0000,0.0000,0.0025),(0.0000,0.0000,0.0000))
Before POSV -- AA=[6,6]((2.0013,0,0,0,0,0),(0,0,0,0,2.0013,0),(0,0,0,2.0013,0,0),(0,0,2.0013,0,0,0),(0,2.0013,0,0,0,0),(2.0013,0,0,0,0,0))
Before POSV -- BB=[6,3]((0.7627,0.1149,0.0025),(0.1149,0.7652,0.1149),(0.0025,0.1149,0.7652),(0,0.0025,0.1149),(0,0,0.0025),(0,0,0))
After POSV -- A=[6,6]((2.0013,0.0000,0.0000,0.0000,0.0000,0.0000),(0.0000,0.0000,0.0000,0.0000,2.0013,0.0000),(0.0000,0.0000,0.0000,2.0013,0.0000,0.0000),(0.0000,0.0000,2.0013,0.0000,0.0000,0.0000),(0.0000,2.0013,0.0000,0.0000,0.0000,0.0000),(2.0013,0.0000,0.0000,0.0000,0.0000,0.0000))
After POSV -- B=[6,3]((0.7627,0.1149,0.0025),(0.1149,0.7652,0.1149),(0.0025,0.1149,0.7652),(0.0000,0.0025,0.1149),(0.0000,0.0000,0.0025),(0.0000,0.0000,0.0000))
After POSV -- AA=[6,6]((1.41467,0,0,0,0,0),(0,0,0,0,2.0013,0),(0,0,0,2.0013,0,0),(0,0,2.0013,0,0,0),(0,2.0013,0,0,0,0),(2.0013,0,0,0,0,0))
After POSV -- BB=[6,3]((0.7627,0.1149,0.0025),(0.1149,0.7652,0.1149),(0.0025,0.1149,0.7652),(0,0.0025,0.1149),(0,0,0.0025),(0,0,0))
So is there another way to use numeric_bindings with plain C array?
Thank you very much.
Best,
-- Marco