February 2012 - petsc-dev - mailman3.cels.anl.gov

Re: [petsc-dev] MatMatMult gives different results
by Hong Zhang 07 Feb '12

07 Feb '12

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, agrayver(a)gfz-potsdam.de < agrayver(a)gfz-potsdam.de> 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" <hzhang(a)mcs.anl.gov> > To: "For users of the development version of PETSc" <petsc-dev(a)mcs.anl.gov > > > 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 < >>>> agrayver(a)gfz-potsdam.de> 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 >> >> >

1 0

Re: [petsc-dev] MatMatMult gives different results
by agrayver＠gfz-potsdam.de 07 Feb '12

07 Feb '12

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" <hzhang(a)mcs.anl.gov> To: "For users of the development version of PETSc" <petsc-dev(a)mcs.anl.gov> 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 <agrayver(a)gfz-potsdam.de> 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

1 0

Re: [petsc-dev] MatMatMult gives different results
by Hong Zhang 07 Feb '12

07 Feb '12

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 < >>> agrayver(a)gfz-potsdam.de> 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 > >

1 0

Re: [petsc-dev] MatMatMult gives different results
by Alexander Grayver 07 Feb '12

07 Feb '12

On 07.02.2012 18:17, Hong Zhang wrote: > Alexander : > >> >> MatCreateMPIAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, >> PETSC_DECIDE, N, M, 12, PETSC_NULL, 12, PETSC_NULL, &AT); >> MatTranspose(A,MAT_IGNORE_MATRIX,&AT); >> >> >> should be replaced with >> MatTranspose(A,MAT_INITIAL_MATRIX,&AT); > Hm, it seems this was the reason... > > With this change, do you still get different matrix product comparing > with Matlab's? > I'll add an error flag in petsc to prevent using > MatTranspose(A,MAT_IGNORE_MATRIX,&AT). Hong, 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. 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! 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 >>> <agrayver(a)gfz-potsdam.de <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

1 0

Re: [petsc-dev] MatMatMult gives different results
by Hong Zhang 07 Feb '12

07 Feb '12

Alexander : > >> MatCreateMPIAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, M, 12, >> PETSC_NULL, 12, PETSC_NULL, &AT); >> MatTranspose(A,MAT_IGNORE_MATRIX,&AT); >> > > should be replaced with > MatTranspose(A,MAT_INITIAL_MATRIX,&AT); > > Hm, it seems this was the reason... > With this change, do you still get different matrix product comparing with Matlab's? I'll add an error flag in petsc to prevent using MatTranspose(A,MAT_IGNORE_MATRIX,&AT). > > 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 < >> agrayver(a)gfz-potsdam.de> 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 > >

1 0

Re: [petsc-dev] MatMatMult gives different results
by Alexander Grayver 07 Feb '12

07 Feb '12

Hong, On 07.02.2012 17:06, Hong Zhang wrote: > Alexander: > Thanks for sending matrix data and testing code. > Before testing it, I have question about your code: > > > MatCreateMPIAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, > M, 12, PETSC_NULL, 12, PETSC_NULL, &AT); > > MatTranspose(A,MAT_IGNORE_MATRIX,&AT); > > > should be replaced with > MatTranspose(A,MAT_INITIAL_MATRIX,&AT); Hm, it seems this was the reason... > > To compare results I use following Matlab script: > > Cm = A.'*B; > > ^ > > 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 .* > > 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 >> <agrayver(a)gfz-potsdam.de <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

1 0

Re: [petsc-dev] MatMatMult gives different results
by Hong Zhang 07 Feb '12

07 Feb '12

Alexander: Thanks for sending matrix data and testing code. Before testing it, I have question about your code: > > MatCreateMPIAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, M, 12, > PETSC_NULL, 12, PETSC_NULL, &AT); > > MatTranspose(A,MAT_IGNORE_MATRIX,&AT); > should be replaced with MatTranspose(A,MAT_INITIAL_MATRIX,&AT); > > To compare results I use following Matlab script: > > Cm = A.'*B; > ^ Do you mean 'Cm = A'*B;'? 'Cm = A.'*B;' gives component-wise matrix product, not matrix product. 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 <agrayver(a)gfz-potsdam.de > > 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 > >

1 0

Re: [petsc-dev] MatMatMult gives different results
by Alexander Grayver 07 Feb '12

07 Feb '12

Dear Hong, You can get matrices here: http://dl.dropbox.com/u/60982984/A.dat http://dl.dropbox.com/u/60982984/B.dat I wrote simple testing application for that problem (it looks ugly, but allows to reproduce the problem): Mat A,AT,B,C; PetscInt M,N; PetscInitialize(&argc,&args,(char *)0,help); ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"A.dat",FILE_MODE_READ,&viewer);CHKERRQ(ierr); ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); ierr = MatSetType(A,MATMPIAIJ);CHKERRQ(ierr); ierr = MatLoad(A,viewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"B.dat",FILE_MODE_READ,&viewer);CHKERRQ(ierr); ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr); ierr = MatSetType(B,MATMPIDENSE);CHKERRQ(ierr); ierr = MatLoad(B,viewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); MatGetSize(A,&M,&N); MatCreateMPIAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, M, 12, PETSC_NULL, 12, PETSC_NULL, &AT); MatTranspose(A,MAT_IGNORE_MATRIX,&AT); MatMatMult(AT,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C); ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"C.dat",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); ierr = PetscViewerSetFormat(viewer,PETSC_VIEWER_NATIVE);CHKERRQ(ierr); ierr = MatView(C,viewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); ierr = PetscFinalize(); return 0; To compare results I use following Matlab script: A = PetscBinaryRead('A.dat','complex',true); B = PetscBinaryRead('B.dat','complex',true); Cm = A.'*B; 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 > <agrayver(a)gfz-potsdam.de <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

1 0

Re: [petsc-dev] [petsc-users] VECMPICUSP with ghosted vector
by Barry Smith 06 Feb '12

06 Feb '12

On Feb 6, 2012, at 3:55 PM, Jed Brown wrote: > On Mon, Feb 6, 2012 at 22:11, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > On the other hand you could argue that ALL vector types should support ghosting as a natural thing (with sequential vectors just have 0 length ghosts conceptually) then it would be desirable to allow setting the ghost information in any ordering. > > Just a comment: Ghost vectors are occasionally meaningful even for serial vectors (if implementing boundary conditions with ghost cells, also coupling). This is fine, currently one needs to use the MPI vector implementation to get this affect; I think it is natural this way and not also require VECSEQ to support this kind of ghosting. Barry

1 0

Re: [petsc-dev] [petsc-users] VECMPICUSP with ghosted vector
by Jed Brown 06 Feb '12

06 Feb '12

On Mon, Feb 6, 2012 at 22:11, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > On the other hand you could argue that ALL vector types should support > ghosting as a natural thing (with sequential vectors just have 0 length > ghosts conceptually) then it would be desirable to allow setting the ghost > information in any ordering. Just a comment: Ghost vectors are occasionally meaningful even for serial vectors (if implementing boundary conditions with ghost cells, also coupling).

1 0