lapackpp-devel

[Lapackpp-devel] LaBandFactDouble Question

[Christophe P. <pi...@tl...>]

Hello,

I am trying to use the banded matrices to solve pde, and I notice that the
return arguments in the function LaBandMatFactorize are weird. I mean that
the pivot and upper/lower matrix should be save in the member B and not in
the original matrix:

So my question is in bfd.h :

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, &A(0,0), &LDA,
&(AF.pivot()(0)), &info);

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.
This is just a question/remark.
Thanks
Christophe

Re: [Lapackpp-devel] LaBandFactDouble Question

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

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

Re: [Lapackpp-devel] LaBandFactDouble Question

[Christophe P. <pi...@tl...>]

Dear Christian,

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 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
>

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
>>
> 
> 

Re: [Lapackpp-devel] LaBandFactDouble Question

[Christophe P. <pi...@tl...>]

Dear Christian,

Thank you for the answer.
Christian Stimming wrote:

> 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.

It is first what I tried, but it seems that the factorize matrix and the 
pivot are not correct ( I compare the results with lapack77 code that is 
proved to be correct). When I print the content of each matrix AF 
doesn't have the correct form, and the pivot vector is completely 
incoherent.
So I tried what you suggested last time, and redesigned a class based on 
the philosophy you used in gfqrc.h. The thing is  I can get the correct 
matrix this time (I mean the correct factorize matrix), but the pivot is 
obviously erroneous (out of bounds references to the matrix). I don't 
think that my design was correct, it was a quick prototype for my 
problem, but the integer pointer to the memory have a strange behavior .
I will keep looking at my implementation of this class, and if I managed 
to make it work, I will be glad to send it to you.
Best regards,
Christophe Picard

> 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