|
|
Ublas : |
Subject: Re: [ublas] [bindings] Possible bug in GBCON
From: Nizar Khalifa Sallem (nksallem_at_[hidden])
Date: 2010-07-27 05:43:44
At Tue, 27 Jul 2010 11:32:28 +0200,
Nizar Khalifa Sallem wrote:
>
> 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
> _______________________________________________
> ublas mailing list
> ublas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: nksallem_at_[hidden]
>
Sorry,
This message was sent by mistake, so I do apologize for all of you :(