Re: [petsc-dev] MatMatMult gives different results

8 Feb 2012

      Hong,

Everything seems to be working well so far.
I would say there is nothing to acknowledge, but I don't mind...

Please let me know if there will be support for mpiaij_mpidense.

Thanks for this great support.

On 08.02.2012 18:43, Hong Zhang wrote:
...
Alexander:
Recently, we add error flags in petsc-dev 
requiring MatXXXSetPreallocation() for many matrix creation routines.
I've updated relevant routines for your code, and add your contributed 
testing code as a petsc example with acknowledge:
petsc-dev/src/mat/examples/tests/ex165.c
Let me know if you do not want your name to be acknowledge in the example.
I tested ex165 with your A and B in sequential (np=1) and parallel 
(np=2, 6):
mpiexec -n <np> ./ex165 -fA A.dat -fB B.dat -view_C
Checking the size of output file C.dat, I see the sequential run and 
parallel run give identical bit size.
Please get the updated petsc-dev and test it. Let me know if you still 
have problem.
It is inefficient to compute C=A^T*B via MatTranspose() and MatMatMult().
I'll see if I can write a C=A^T*B  for mpiaij_mpidense.
Hong
In addition, I changed my code and form A^T explicitly (to avoid
    previous error) and I got another one:
[1]PETSC ERROR: --------------------- Error Message
    ------------------------------------
    [1]PETSC ERROR: Object is in wrong state!
    [1]PETSC ERROR: Must call MatXXXSetPreallocation() or MatSetUp()
    on argument 1 "mat" before MatAssemblyBegin()!
    [1]PETSC ERROR:
    ------------------------------------------------------------------------
    [1]PETSC ERROR: Petsc Development HG revision:
    249597282bcb6a1051042a9fdfa5705679ed4f18  HG Date: Tue Feb 07
    09:44:23 2012 -0600
    [1]PETSC ERROR: See docs/changes/index.html for recent updates.
    [1]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
    [1]PETSC ERROR: See docs/index.html for manual pages.
    [1]PETSC ERROR:
    ------------------------------------------------------------------------
    [1]PETSC ERROR: solveTest on a openmpi-i named glic1 by agrayver
    Wed Feb  8 13:06:45 2012
    [1]PETSC ERROR: Libraries linked from
    /home/lib/petsc-dev/openmpi-intel-complex-debug-f-mkl/lib
    [1]PETSC ERROR: Configure run at Tue Feb  7 18:19:58 2012
    [1]PETSC ERROR: Configure options
    --with-petsc-arch=openmpi-intel-complex-debug-f-mkl
    --with-fortran-interfaces=1 --download-superlu
    --download-superlu_dist --download-mumps --download-parmetis
    --download-ptscotch --download-metis
    --with-scalapack-lib=/opt/intel/Compiler/11.1/072/mkl/lib/em64t/libmkl_scalapack_lp64.a
    --with-scalapack-include=/opt/intel/Compiler/11.1/072/mkl/include
    --with-blacs-lib=/opt/intel/Compiler/11.1/072/mkl/lib/em64t/libmkl_blacs_openmpi_lp64.a
    --with-blacs-include=/opt/intel/Compiler/11.1/072/mkl/include
    --with-mpi-dir=/opt/mpi/intel/openmpi-1.4.2
    --with-scalar-type=complex
    --with-blas-lapack-lib="[/opt/intel/Compiler/11.1/072/mkl/lib/em64t/libmkl_intel_lp64.a,/opt/intel/Compiler/11.1/072/mkl/lib/em64t/libmkl_intel_thread.a,/opt/intel/Compiler/11.1/072/mkl/lib/em64t/libmkl_core.a,/opt/intel/Compiler/11.1/072/lib/intel64/libiomp5.a]"
    --with-precision=double --with-x=0
    [1]PETSC ERROR:
    ------------------------------------------------------------------------
    [1]PETSC ERROR: MatAssemblyBegin() line 4795 in
    /home/lib/petsc-dev/src/mat/interface/matrix.c
    [1]PETSC ERROR: MatMatMultSymbolic_MPIAIJ_MPIDense() line 638 in
    /home/lib/petsc-dev/src/mat/impls/aij/mpi/mpimatmatmult.c
    [1]PETSC ERROR: MatMatMult_MPIAIJ_MPIDense() line 594 in
    /home/lib/petsc-dev/src/mat/impls/aij/mpi/mpimatmatmult.c
    [1]PETSC ERROR: MatMatMult() line 8618 in
    /home/lib/petsc-dev/src/mat/interface/matrix.c
On 08.02.2012 17:45, Hong Zhang wrote:
...
Alexander :
    I can repeat the crash, and am working on it.
    I'll let you know after the bug is fixed.
Thanks for your patience,
    Hong
It seems now I know why I used MAT_REUSE_MATRIX and
        preallocated matrix.
        I removed MatCreateMPIAIJ() and use just MatTranspose with
        MAT_INITIAL_MATRIX.
        As a consequence I get hundreds of errors like:
[30]PETSC ERROR: MatSetValues_MPIAIJ() line 538 in
        /home/mt/agrayver/lib/petsc-dev/src/mat/impls/aij/mpi/mpiaij.c
        [30]PETSC ERROR: MatAssemblyEnd_MPIAIJ() line 653 in
        /home/mt/agrayver/lib/petsc-dev/src/mat/impls/aij/mpi/mpiaij.c
        [30]PETSC ERROR: MatAssemblyEnd() line 4978 in
        /home/mt/agrayver/lib/petsc-dev/src/mat/interface/matrix.c
        [30]PETSC ERROR: MatTranspose_MPIAIJ() line 2061 in
        /home/mt/agrayver/lib/petsc-dev/src/mat/impls/aij/mpi/mpiaij.c
        [30]PETSC ERROR: MatTranspose() line 4397 in
        /home/mt/agrayver/lib/petsc-dev/src/mat/interface/matrix.c
        [31]PETSC ERROR: --------------------- Error Message
        ------------------------------------
        [31]PETSC ERROR: Argument out of range!
        [31]PETSC ERROR: New nonzero at (1659,53337) caused a malloc!
Ans this is the lastest petsc-dev revision.
        I know there were some changes in petsc-dev concerning this
        issue.
        Can you give me a hint how to avoid this?
As for why I need C=A^T*B. This product is used further as a
        system matrix for LSQR solver. The largest dimension is on
        the order of 10^6
        Since I form A and B myself I can, of course, form A^T
        explicitly, but I thought I would first implement everything
        as it is written down on paper and then optimize once it
        works (I thought it's easier way).
On 07.02.2012 20:44, Hong Zhang wrote:
...
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]
        <mailto:[email protected]> <[email protected]
        <mailto:[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]
            <mailto:[email protected]>>
            To: "For users of the development version of PETSc"
            <[email protected] <mailto:[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]
                        <mailto:[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
-- 
        Regards,
        Alexander
-- 
    Regards,
    Alexander
-- 
Regards,
Alexander

Alexander Grayver

tags

participants (1)