Alexander, I'm curious about why do you need parallel C=A^T*B? How large your matrices are? In petsc-dev, we have MatTransposeMatMult() for mpiaij and mpiaij, but not mpiaij and mpidense. We may add support of MatTransposeMatMult_MPIAIJ_MPIDense() if there is such need. Hong On Tue, Feb 7, 2012 at 1:18 PM, [email protected] < [email protected]> wrote:

Hong,

Thanks for explanation. I will try this tomorrow. Good to have this stuff in the help now.

And sorry for misleading you initially.

Regards, Alexander

----- Reply message ----- From: "Hong Zhang" <[email protected]> To: "For users of the development version of PETSc" <[email protected]

...

Subject: [petsc-dev] MatMatMult gives different results Date: Tue, Feb 7, 2012 19:09

Alexander :

...

There is something I didn't get yet, I hope you could clarify it.

So, when I use flag MAT_INITIAL_MATRIX in test program it works fine.

Good to know :-)

...

If I put this flag in my original program I get dozens of exceptions like: [42]PETSC ERROR: Argument out of range! [42]PETSC ERROR: New nonzero at (1336,153341) caused a malloc!

You cannot do MatCreateMPIAIJ() MatTranspose(A,MAT_INITIAL_MATRIX,&AT);

MatCreateMPIAIJ() creates AT and preallocates approximate nonzeros, which does not match exactly the nonzeros in MatTranspose(A,MAT_INITIAL_MATRIX,&AT); MatTranspose(A,MAT_INITIAL_MATRIX,&AT) creates matrix AT and sets correct nonzero pattern and values in AT. MatTranspose() only takes in "MAT_INITIAL_MATRIX" - for a new AT, and "MAT_REUSE_MATRIX" when AT is created with MatTranspose(A,MAT_INITIAL_MATRIX,&AT) and reuse for updating its values (not nonzero patten).

I'm updating petsc help menu on MatTranspose(). Thanks for the report.

Hong

...

I changed this flag to MAT_REUSE_MATRIX and exceptions disappeared, but result is incorrect again (same as for MAT_IGNORE_MATRIX) I tried test program with MAT_REUSE_MATRIX and it also gives different matrix product.

Since there is no description of MatReuse structure for MatTranspose it's a bit confusing what to expect from it.

...

Do you mean 'Cm = A'*B;'? 'Cm = A.'*B;' gives component-wise matrix product, not matrix product.

.' operator means non-Hermitian transpose. That is what I get with MatTranspose (in contrast with MatHermitianTranspose) component-wise matrix product would be .*

You are correct.

Hong

...

Hong

...

C = PetscBinaryRead('C.dat','complex',true);

Matrix C is different depending on number of cores I use. My PETSc is: Using Petsc Development HG revision: 876c894d95f4fa6561d0a91310ca914592527960 HG Date: Tue Jan 10 19:27:14 2012 +0100

On 06.02.2012 17:13, Hong Zhang wrote:

MatMatMult() in petsc is not well-tested for complex - could be buggy. Can you send us the matrices A and B in petsc binary format for investigation?

Hong

On Mon, Feb 6, 2012 at 5:55 AM, Alexander Grayver < [email protected]> wrote:

...

Dear PETSc team,

I try to use: call MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT_DOUBLE_PRECISION,C,ierr);CHKERRQ(ierr)

Where both A and B are rectangular, but A is sparse and B is dense. Both are double complex and distributed. The product PETSc gives me contains some errors in some part of the matrix. I output A, B and C then computed product in matlab.

Attached you see figure plotted as: imagesc(log10(abs(C-Cm)))

Where Cm -- product computed in matlab.

The pattern and amplitude vary depending on the number of cores I use. This picture is obtained for 48 cores (I've tried 12, 64 cores as well).

Where should I look for possible explanation?

-- Regards, Alexander

-- Regards, Alexander

-- Regards, Alexander

-- Regards, Alexander