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 >
