Question on MatPtAP

[Hong Zhang]

MatPtAP() should work for rectangular P.
I'll test your code and let you know what I get,
hopefully later this evening,

Hong

On Tue, 1 Jul 2008, Tobias Neckel wrote:

> Hello,
>
> when moving from PETSc 2.3.2-p10 to 2.3.3-p13, I recently encountered a 
> strange behaviour (segmentation fault) of the new version concerning the 
> usage of MatPtAP. I tried to pull everything down to a simple test case. For 
> square matrices (A as well as P), everything seems to work fine.
>
> But when I set up a simple non-square example (A being the 3x3-identity, P 
> being a 3x5 matrix, see attached file PETScLibTest.cpp), I encounter severe 
> problems: As soon as P gets entries outside its 3x3 block (an entry P(0,3), 
> e.g.), PETSc is telling me about Memory corruption while using MatPtAP (see 
> part of command line output in the attached file output.txt, run with -info).
>
> I used valgrind to check if sth. weird happens, but it is giving nothing 
> until the PETSc error message.
>
> Using the same test case with PETSc 2.3.2-p10 (with identic configuring 
> options on the same machine) does not show any problems at all.
>
> So I am a bit confused. Is it not allowed any longer to use MatPtAP with 
> non-square matrices? I checked the online documentation but did not find 
> anything ...
>
> Next idea was to use MatMatMult directly to see if that works. So I used
> MatMatMultTranspose(P, A, MAT_INITIAL_MATRIX, 1.0, &C1) and then
> MatMatMult(C1, P, MAT_INITIAL_MATRIX, 1.0, &C)
> to get C=C1*P=P^T*A*P in two steps. This time, everything works fine for the 
> small test case from above, also with the new version 2.3.3-p13.
>
>
> Best regards
> Tobias Neckel
>
> -- 
> Dipl.-Tech. Math. Tobias Neckel
>
> Institut für Informatik V, TU München
> Boltzmannstr. 3, 85748 Garching
>
> Tel.:   089/289-18602
> Email:  neckel at in.tum.de
> URL:    http://www5.in.tum.de/persons/neckel.html
>