|
|
Ublas : |
Subject: Re: [ublas] [bindings] Possible bug in GBCON
From: Nizar Khalifa Sallem (nksallem_at_[hidden])
Date: 2010-07-27 05:32:28
At Tue, 27 Jul 2010 11:00:50 +0200,
Marco Guazzone wrote:
>
> [1 <text/plain; ISO-8859-1 (7bit)>]
> Hi,
>
> I'm facing a problem with GBCON.
>
> I use it to estimate the reciprocal matrix condition number of a banded matrix.
>
> Firstly, I call GBTRF to get the LUP decomposition (AB, IPIV) of my
> input matrix A.
> Then, I compute the matrix norm (specifically, the 1-norm).
> Finally, I call GBCON using as parameters the above (AB,IPIV) pair and
> the matrix norm:
> bindings::lapack::gbcon('O', AB, IPIV, norm, rcond)
> The result is stored in 'rcond'.
>
> Well, when I call GBCON the following assertion fails (gbcon.hpp -
> lines 155,250):
>
> BOOST_ASSERT( bindings::stride_major(ab) >=
>
> 2*bindings::bandwidth_lower(ab)+bindings::bandwidth_upper(ab)+1 );
>
> I think the problem is in the use of
> bindings::bandwidth_upper(ab)
> which IMHO should be replaced with
> bindings::bandwidth_upper(ab) - bindings::bandwidth_lower(ab)
> As a matter of fact, LAPACK doc
> http://www.netlib.org/lapack/explore-html/a00579_source.html
> refers to KL and KU with respect to the original input matrix A and
> not to its LUP decomposition AB.
>
> I submit a test case and a possible patch (the patch also fixes a
> similar issue when calling detail::gbcon).
>
> What do you think?
>
> Thank you very much!!
>
> Cheers,
>
> -- Marco
> [2 lapack_gbcon.cpp <text/x-c++src; US-ASCII (base64)>]
> #include <algorithm>
> #include <boost/numeric/bindings/lapack/computational/gbcon.hpp>
> #include <boost/numeric/bindings/lapack/computational/gbtrf.hpp>
> #include <boost/numeric/bindings/lapack/auxiliary/lange.hpp>
> #include <boost/numeric/bindings/ublas.hpp>
> #include <boost/numeric/ublas/banded.hpp>
> #include <boost/numeric/ublas/io.hpp>
> #include <boost/numeric/ublas/matrix.hpp>
> #include <boost/numeric/ublas/matrix_expression.hpp>
> #include <boost/numeric/ublas/operation/num_columns.hpp>
> #include <boost/numeric/ublas/operation/num_rows.hpp>
> #include <boost/numeric/ublas/traits.hpp>
> #include <boost/numeric/ublas/vector.hpp>
> #include <cstddef>
> #include <iostream>
> #include <string>
>
>
> namespace ublas = boost::numeric::ublas;
> namespace bindings = boost::numeric::bindings;
>
>
> int main()
> {
> typedef double value_type;
> typedef ublas::banded_matrix<value_type, ublas::column_major> matrix_type;
> typedef ublas::matrix_traits<matrix_type>::size_type size_type;
> typedef matrix_type::array_type array_type;
> typedef ublas::type_traits<value_type>::real_type result_type;
> typedef ublas::banded_matrix<value_type,ublas::row_major,array_type> auxiliary_matrix_type;
>
> const size_type nr = 4;
> const size_type nc = 4;
> const size_type kl = 1;
> const size_type ku = 2;
> const size_type n = std::min(nr, nc);
>
> matrix_type A(nr,nc,kl,ku);
> A(0,0) = -0.23; A(0,1) = 2.54; A(0,2) = -3.66; /* 0 */ // 1st row
> A(1,0) = -6.98; A(1,1) = 2.46; A(1,2) = -2.73; A(1,3) = -2.13; // 2nd row
> /* 0 */ A(2,1) = 2.56; A(2,2) = 2.46; A(2,3) = 4.07; // 3rd row
> /* 0 */ /* 0 */ A(3,2) = -4.78; A(3,3) = -3.82; // 4th row
>
>
> // Compute LUP decomposition
>
> auxiliary_matrix_type AB(A, kl, kl+ku);
>
> std::cout << "AB = " << AB << std::endl;
>
> ublas::vector< ::fortran_int_t > ipiv(n);
>
> bindings::lapack::gbtrf(AB, ipiv);
>
> std::cout << "GBTRF -> AB = " << AB << std::endl;
> std::cout << "GBTRF -> IPIV = " << ipiv << std::endl;
>
> // Compute matrix 1-norm
>
> result_type norm(0);
>
> //NOTE: bindings::lapack::lange is actually broken
>
> norm = ublas::norm_1(A);
>
> // Estimate reciprocal matrix condition number
>
> result_type rcond(0);
>
> bindings::lapack::gbcon('O', AB, ipiv, norm, rcond);
>
> std::cout << "GBCON -> RCOND = " << rcond << std::endl;
> }
> [3 gbcon.patch <application/octet-stream (base64)>]
>
> [4 <text/plain; us-ascii (7bit)>]
> _______________________________________________
> ublas mailing list
> ublas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: nksallem_at_[hidden]
Salut,
J'ai identifié deux problèmes le premier c'est boost qui n'est pas a
la bonne version ensuite tu devras faire un svn up dans bundler parce
que ton CMakeLists est un quelque peu erroné, c'est de ma faute je ne
l'ai pas fait. Je te demande pardon.
Pour boost je vais voir si tu as un robotpkg quelque part sinon je
vais te l'installer ou voir comment tu as fait