Pierre, I believe you are in the correct direction for debugging MatProductReplaceMats() . I 'll investigate it and let you know the result. Hong
________________________________ From: Pierre Jolivet <pie...@joliv.et> Sent: Friday, July 15, 2022 12:01 AM To: Zhang, Hong <hzh...@mcs.anl.gov> Cc: Barry Smith <bsm...@petsc.dev>; For users of the development version of PETSc <petsc-dev@mcs.anl.gov> Subject: Re: [petsc-dev] MatProduct_AtB --with-scalar-type=complex Barry, MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ() is indeed called. product->alg is default, square is PETSC_FALSE. Hong, I believe the issue comes from the fact that atb->updateAt is PETSC_FALSE in MatProductNumeric_AtB_SeqAIJ_SeqAIJ(). If the name of this variable is relevant to its purpose, I believe it should be set to PETSC_TRUE when calling MatProductReplaceMats() whenever A is changed. I would prefer using MatProductReplaceMats() because I’m implementing the same MatConvert() as MatNormal for the Hermitian case and it’s the only way to reuse the symbolic product, cf. https://petsc.org/main/src/mat/impls/normal/normm.c.html#line315 in the case where --with-scalar-type=real Thanks, Pierre On 15 Jul 2022, at 5:52 AM, Zhang, Hong <hzh...@mcs.anl.gov<mailto:hzh...@mcs.anl.gov>> wrote: Pierre, Our MatProductReplaceMats() is not well tested, which might be buggy. I simplified your code without calling MatProductReplaceMats() and got correct results in the cases ./ex1111 -product_view ::ascii_matlab -convert false/true -correct false and ./ex1111 -product_view ::ascii_matlab -convert false/true -correct true My code is attached. I'll investigate MatProductReplaceMats(). Hong ________________________________ From: petsc-dev <petsc-dev-boun...@mcs.anl.gov<mailto:petsc-dev-boun...@mcs.anl.gov>> on behalf of Barry Smith <bsm...@petsc.dev<mailto:bsm...@petsc.dev>> Sent: Thursday, July 14, 2022 4:38 PM To: Pierre Jolivet <pie...@joliv.et<mailto:pie...@joliv.et>> Cc: For users of the development version of PETSc <petsc-dev@mcs.anl.gov<mailto:petsc-dev@mcs.anl.gov>> Subject: Re: [petsc-dev] MatProduct_AtB --with-scalar-type=complex Can you confirm if MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ() ends up being called for you and what path it takes inside that routine (depends) on the algorithm it is using. > On Jul 14, 2022, at 4:30 PM, Pierre Jolivet > <pie...@joliv.et<mailto:pie...@joliv.et>> wrote: > > Hello, > In the following example, the SeqAIJ implementation of MatProduct_AtB produce > a different (and wrong) result, compared to the SeqDense implementation or > MATLAB. > I want to compute B = A^H A (where ^H is the Hermitian transpose). > So I create a MatProduct with A and A. > Duplicate A into another Mat which I conjugate. > And I replace the first Mat of the product with this conjugate. > I expect to get the proper result, which I don’t. > Is the MatProduct_AtB implementation in the complex case not computing A^T B > (where ^T is the transpose)? > For reference, here is how to properly compute A^H A with current main: > conj(A^H conj(A)) — so it requires an extra MatConjugate I’d like to avoid. > > Thanks, > Pierre > > <ex1111.c> > > $ ./ex1111 -product_view ::ascii_matlab -A_view ::ascii_matlab -convert false > %Mat Object: 1 MPI process > % type: seqdense > % Size = 2 2 > Mat_0xc4000001_0 = zeros(2,2); > Mat_0xc4000001_0 = [ > 7.2003197397953400e-01 + 6.1793966542126100e-02i 3.9777780919128602e-01 + > 7.3036588248200474e-02i > 1.0022337819588500e-02 + 1.4463931936456476e-01i 1.0386628927366459e-01 + > 2.5078039364333193e-01i > ]; > %Mat Object: 1 MPI process > % type: seqdense > % Size = 2 2 > Mat_0xc4000001_1 = zeros(2,2); > Mat_0xc4000001_1 = [ > 5.4328551781548817e-01 + 0.0000000000000000e+00i 3.2823965013353340e-01 + > 1.5498666614872689e-02i > 3.2823965013353340e-01 + -1.5498666614872689e-02i 2.3724054059134142e-01 + > 0.0000000000000000e+00i > ]; > > $ ./ex1111 -product_view ::ascii_matlab -convert true > %Mat Object: 1 MPI process > % type: seqaij > % Size = 2 2 > % Nonzeros = 4 > zzz = zeros(4,4); > zzz = [ > 1 1 4.9380746380098023e-01 9.1886511660038694e-02 > 1 2 2.4666779825931440e-01 9.4705502650537468e-02 > 2 1 2.4666779825931440e-01 9.4705502650537468e-02 > 2 2 1.0079024247365802e-01 1.1019992594899400e-01 > ]; > Mat_0xc4000001_0 = spconvert(zzz); > > $ ./ex1111 -product_view ::ascii_matlab -convert true -correct true > %Mat Object: 1 MPI process > % type: seqaij > % Size = 2 2 > % Nonzeros = 4 > zzz = zeros(4,4); > zzz = [ > 1 1 5.4328551781548828e-01 -0.0000000000000000e+00 > 1 2 3.2823965013353340e-01 1.5498666614872696e-02 > 2 1 3.2823965013353340e-01 -1.5498666614872696e-02 > 2 2 2.3724054059134142e-01 -0.0000000000000000e+00 > ]; > Mat_0xc4000001_0 = spconvert(zzz); > > <Screenshot 2022-07-14 at 10.12.53 PM.png> <ex1111.c>
