Re: [ublas] Fwd: boost/numeric/bindings/lapack

Subject: Re: [ublas] Fwd: boost/numeric/bindings/lapack
From: Vardan Akopian (vakopian_at_[hidden])
Date: 2008-11-13 20:49:58

On Thu, Nov 13, 2008 at 1:13 AM, Julius Ziegler <ziegler_at_[hidden]> wrote:
> Andreas Klöckner wrote:
>>
>> On Mittwoch 12 November 2008, Thomas Klimpel wrote:
>>
>>>>
>>>> I've received the following patch suggestion from Julius Ziegler. What
>>>> he
>>>> points out does look fishy. Since I've never worked with banded matrices
>>>> in LAPACK, I hope someone on the list is better equipped to give a
>>>> definitive answer.
>>>>
>>>
>>> The fix from Julius Ziegler is correct (according to the documentation of
>>> SGBTRF and ublas/banded_matrix), but it doesn't go far enough. It should
>>> also change "assert(ku > 0)" into "assert(ku >= 0)", because that is all
>>> that is required by SGBTRF.
>>>
>>
>> I've attached a patch. Let me know if it looks ok, and if so, I'll apply
>> it to my tree and release a new snapshot. Note that I've also removed the
>> comments above the code snippets in question--they don't seem to apply any
>> more. Tell me if I'm wrong.
>>
>> Andreas
>>
>
> I believe now that my suggestion is wrong. SGBTRF seems to require kl extra
> diagonals allocated within the matrix band, so the code is correct.
> ublas::sgbrtf does not give correct results for me, though, I will try to
> put together an example an post it here.
>
> Cheers,
> Julius

Julius,
Please make sure that you don't have BOOST_UBLAS_OWN_BANDED defined.
Also make sure that you use row_major as the layout for the
banded_matrix (this is the default, and is also checked within
gbtrf(), if BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK is not defined).
Also make sure that you actually do allocate kl+ku upper diagonals in
your banded matrix. So if you want a 100x100 banded matrix with 3
lower and 2 upper diagonals, you should do:
ublas::banded_matrix<double> a(100, 100, 3, 2+3);
// fill the values of a ...
ublas::vector<int> piv(100);
int info = lapack::gbtrf(a, piv);
if (info != 0) {
    // error
}

-Vardan