Re: [Lapackpp-devel] LaBandFactDouble Question

[Christian S. <sti...@tu...>]

Dear Christophe,

Christophe Picard schrieb:
> Thank you for your answer.
> Actually, I am not really a good programmer, but I am quite familiar 
> with the usage of matrices for solving PDE. Actually, a lot of people 
> who are using LAPACK F77 use the banded matrix vs the full matrix for 
> the memory efficiency. Personally, when I use it, I overwrite the 
> original matrix since the only reason to use LU is because the solving 
> part is really fast when compare to GMRES or other iterative methods, 
> and the overhead of computing the factorization it is becoming 
> negligible when the solver is used a lot.

I totally agree that this solving method is very efficient in specific 
situations, and I would certainly like to offer the correct program code 
in lapackpp.

After looking at this in more detail, I'd say that what you wanted to do 
can be achieved by constructing the LaBandFactDouble object with the 
original matrix A as argument, so that the original matrix A is used as 
storage for the decomposition/factorization?

LaBandMatDouble A(m,n);
// fill A ..., then:
LaBandFactDouble AF(A);
LaBandMatFactorize(A, AF);
// A and AF both refer to the factorization, which is fine:
LaLinearSolve(AF, X, B); // ... and so on.

This does not mean that the current lapackpp code is correct in any 
case, but at least it should work for you. The point here is that the 
original matrix A is no longer being used after constructing the 
factorization.

The only real solution would be a full redesign of the factorization 
classes, somehow along the lines of the gfqrc.h. But for now, the above 
code should work for you.

Best regards,

Christian Stimming


> I will try to see if I can write something. I am pretty interested in 
> solver, and I would really like to make this one working.
> 
> Regards,
> Christophe
> 
> 
> Christian Stimming wrote:
> 
>> Dear Christophe,
>>
>> Christophe Picard schrieb:
>>
>>> I am trying to use the banded matrices to solve pde, and I notice 
>>> that the
>>> return arguments in the function LaBandMatFactorize are weird. 
>>
>>
>>
>> Absolutely. Yes, in fact *all* the *FactDouble classes are pretty much 
>> messed up. In gfd.h, I added a comment which made this clear, but in 
>> the other Fact headers I forgot to add such a comment (will be 
>> included immediately).
>>
>> In other words, the design of these classes is screwed anyway. If you 
>> found a way to make them useable, then you would be very welcome to 
>> contribute your code to lapackpp. I will happily throw out the wrong 
>> code and replace it by your code.
>>
>>> I mean that
>>> the pivot and upper/lower matrix should be save in the member B and 
>>> not in
>>> the original matrix:
>>
>>
>>
>> Err, yes, but the whole class needs more thought. It depends on how 
>> you construct the LaBandFactDouble object. If you use the constructor 
>> LaBandFactDouble(LaBanMatDouble& A), then the member variable B 
>> references the same memory as A and the function LaBandMatFactorize 
>> will use the matrix A in-place, regardless whether A or B is passed as 
>> parameter.
>>
>> I'm not sure how such a factorization is used in practical numerics. 
>> It might well be possible that you should throw away the whole 
>> implementation of bfd.h and invent your own. On the other hand, I 
>> recently added the class LaGenQRFactComplex in gfqrc.h and I think its 
>> design is one possibility to efficiently model matrix factorization as 
>> a class. Maybe you can try to copy the implementation model of gfqrc.h 
>> into one that uses LaBandMatDouble, but that's just one possibility 
>> out of many.
>>
>>> So my question is in bfd.h : (...)
>>> should it be replace with something more along this (even maybe overload
>>> the function)
>>>
>>> 00138 inline void LaBandMatFactorize(LaBandMatDouble &A, 
>>> LaBandFactDouble
>>> &AF)
>>> 00139 {
>>> 00140     integer n = A.size(1), m = n, LDA = A.gdim(0);
>>> 00141     integer KL = A.subdiags(), KU = A.superdiags(), info=0;
>>> 00142
>>> 00143     F77NAME(dgbtrf)(&m, &n, &KL, &KU, &AF.B(0,0), &LDA,
>>> &(AF.pivot()(0)), &info);
>>>
>>> otherwise the factorization is not saved into AF but in A, making AF.B
>>> useless, and the function LaLinearSolve unusable.
>>
>>
>>
>> Absolutely. Throw out the old errornous function and replace it by 
>> your own function, and I would happily include that in the official 
>> lapackpp repository. Especially if you have tested and verified your 
>> implementation.
>>
>>> This is just a question/remark.
>>
>>
>>
>> Thank you very much for this feedback.
>>
>> Regards,
>>
>> Christian Stimming
>>
> 
>