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> on behalf of Barry Smith
<bsm...@petsc.dev>
Sent: Thursday, July 14, 2022 4:38 PM
To: Pierre Jolivet <pie...@joliv.et>
Cc: For users of the development version of PETSc <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> 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>
#include <petsc.h>
static char help[] = "";
int main(int argc,char **args)
{
PetscInt n = 2;
Mat array[1],B,conjugate;
PetscBool flg = PETSC_FALSE, correct = PETSC_FALSE;
PetscCall(PetscInitialize(&argc,&args,NULL,help));
PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,n,n,n,n,-1.0,array));
PetscCall(MatConvert(array[0],MATDENSE,MAT_INPLACE_MATRIX,array));
PetscCall(MatSetRandom(array[0],NULL));
PetscCall(MatViewFromOptions(array[0],NULL,"-A_view"));
PetscCall(PetscOptionsGetBool(NULL,NULL,"-convert",&flg,NULL));
PetscCall(PetscOptionsGetBool(NULL,NULL,"-correct",&correct,NULL));
if (flg) {
PetscCall(MatConvert(array[0],MATAIJ,MAT_INPLACE_MATRIX,array));
//PetscCall(PetscOptionsGetBool(NULL,NULL,"-correct",&correct,NULL));
}
PetscCall(MatDuplicate(array[0], MAT_COPY_VALUES, &conjugate));
PetscCall(MatConjugate(conjugate));
if (!correct) PetscCall(MatProductCreate(array[0],array[0],NULL,&B));
else PetscCall(MatProductCreate(conjugate,array[0],NULL,&B));
PetscCall(MatProductSetType(B,MATPRODUCT_AtB));
PetscCall(MatProductSetFromOptions(B));
PetscCall(MatProductSymbolic(B));
PetscCall(MatProductNumeric(B));
PetscCall(MatDestroy(&conjugate));
PetscCall(MatViewFromOptions(B,NULL,"-product_view"));
PetscCall(MatDestroy(&B));
PetscCall(MatDestroy(array));
PetscCall(PetscFinalize());
return 0;
}
