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

Re: [petsc-dev] Performance of VecMDot_SeqCUSP
by Daniel Lowell 24 Apr '12

24 Apr '12

I'm writing a vector type with uses flag synching like you have in PETSc with Vec CUSP, however it uses asynchronous kernel launches (pipeling,etc..) and autotuned kernels. Not quite ready for primetime, but we have seen the value of it in terms of speed up. On Tue, Apr 24, 2012 at 2:09 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > On Tue, Apr 24, 2012 at 11:55, Jose E. Roman <jroman(a)dsic.upv.es> wrote: > >> It seems that VecMDot_SeqCUSP has rather poor performance. This has a lot >> of impact in SLEPc because it is the main kernel used in the >> orthogonalization of vectors. >> >> Is this due to the version of Thrust? I am using CUDA Toolkit 4.0. >> >> I tried a naive replacement that copies the contents of the vectors into >> a matrix and calls CUBLAS dgemv. The improvement is significant, despite >> the data movement overhead. In some tests I see a reduction of time >> (VecReduceArith) from 24.5 seconds to 9.6 seconds (with up to 200 vectors >> of length 10000) on a Fermi. >> >> I can send the code for you to try. >> > > That would be useful. The current code is horribly over-synchronous since > it spawns a new kernel for each set of 3 vectors. We want to use one kernel > to process all vectors. Unless there is another implementation somewhere, > we'll probably want to write custom CUDA kernel for this (to avoid needing > to copy all the vectors into contiguous memory). > -- Daniel Lowell

1 0

Re: [petsc-dev] Performance of VecMDot_SeqCUSP
by Jed Brown 24 Apr '12

24 Apr '12

On Tue, Apr 24, 2012 at 11:55, Jose E. Roman <jroman(a)dsic.upv.es> wrote: > It seems that VecMDot_SeqCUSP has rather poor performance. This has a lot > of impact in SLEPc because it is the main kernel used in the > orthogonalization of vectors. > > Is this due to the version of Thrust? I am using CUDA Toolkit 4.0. > > I tried a naive replacement that copies the contents of the vectors into a > matrix and calls CUBLAS dgemv. The improvement is significant, despite the > data movement overhead. In some tests I see a reduction of time > (VecReduceArith) from 24.5 seconds to 9.6 seconds (with up to 200 vectors > of length 10000) on a Fermi. > > I can send the code for you to try. > That would be useful. The current code is horribly over-synchronous since it spawns a new kernel for each set of 3 vectors. We want to use one kernel to process all vectors. Unless there is another implementation somewhere, we'll probably want to write custom CUDA kernel for this (to avoid needing to copy all the vectors into contiguous memory).

1 0

Re: [petsc-dev] MetSetBlockSize issue
by Mark F. Adams 24 Apr '12

24 Apr '12

On Apr 24, 2012, at 2:55 PM, Barry Smith wrote: > > On Apr 24, 2012, at 1:39 PM, Mark F. Adams wrote: > >> Now I'm getting this error with code like this: >> >> ierr = MatGetSubMatrix( Cmat, new_eq_indices, new_eq_indices, MAT_INITIAL_MATRIX, &mat ); >> CHKERRQ(ierr); >> ierr = MatSetBlockSize( mat, cbs ); CHKERRQ(ierr); >> >> and like this: >> >> ierr = MatPtAP( Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, &Cmat ); CHKERRQ(ierr); >> ierr = MatSetBlockSize( Cmat, cbs ); CHKERRQ(ierr); > > Right you cannot do this. The matrix already exists so you cannot now set its size. > > Does Cmat have a block size of cbs > > What about Amat? By default these routines should create the new matrix with the correct blocksize. > > The harder part is if the blocksize of mat would be different than Cmat? > It seems like MatGetSubMatrix should just inherit the block size but PtAP is harder. The block size of the Cmat is the column block size of P. But I can not set a column block size (!= row block size) so I don't set block size on P at all. Here cbs is this column block size of P, or what it should be. Mark > Barry > >> >> Mark >> >> On Apr 23, 2012, at 9:09 PM, Barry Smith wrote: >> >>> >>> Yes, look, for example how ex32 is handled at the bottom of the makefile. >>> >>> Thanks >>> >>> Barry >>> >>> On Apr 23, 2012, at 5:50 PM, Mark F. Adams wrote: >>> >>>> >>>> On Apr 23, 2012, at 5:59 PM, Barry Smith wrote: >>>> >>>>> >>>>> It was changed a while ago that MatSetBlockSize() couldn't be set after the matrix was full instantiated. I guess that example did not get fixed because it is not listed in the makefile to run in the makeall. >>>> >>>> May I add it? >>>> >>>>> >>>>> I have fixed the example to call MatSetBlockSize() at a safe point. >>>>> >>>>> Barry >>>>> >>>>> On Apr 23, 2012, at 4:29 PM, Mark F. Adams wrote: >>>>> >>>>>> ex55.c in ksp is failing with: >>>>>> >>>>>> [0]PETSC ERROR: --------------------- Error Message ------------------------------------ >>>>>> [0]PETSC ERROR: Arguments are incompatible! >>>>>> [0]PETSC ERROR: Cannot change block size 1 to 2! >>>>>> [0]PETSC ERROR: ------------------------------------------------------------------------ >>>>>> >>>>>> on this line 57 of ex55.c: >>>>>> >>>>>> ierr = MatSetBlockSize(Amat,2); CHKERRQ(ierr); >>>>>> >>>>>> Any idea what happened here? >>>>>> >>>>>> Mark >>>>> >>>>> >>>> >>> >>> >> > >

1 0

Re: [petsc-dev] MetSetBlockSize issue
by Barry Smith 24 Apr '12

24 Apr '12

On Apr 24, 2012, at 1:39 PM, Mark F. Adams wrote: > Now I'm getting this error with code like this: > > ierr = MatGetSubMatrix( Cmat, new_eq_indices, new_eq_indices, MAT_INITIAL_MATRIX, &mat ); > CHKERRQ(ierr); > ierr = MatSetBlockSize( mat, cbs ); CHKERRQ(ierr); > > and like this: > > ierr = MatPtAP( Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, &Cmat ); CHKERRQ(ierr); > ierr = MatSetBlockSize( Cmat, cbs ); CHKERRQ(ierr); Right you cannot do this. The matrix already exists so you cannot now set its size. Does Cmat have a block size of cbs What about Amat? By default these routines should create the new matrix with the correct blocksize. The harder part is if the blocksize of mat would be different than Cmat? Barry > > Mark > > On Apr 23, 2012, at 9:09 PM, Barry Smith wrote: > >> >> Yes, look, for example how ex32 is handled at the bottom of the makefile. >> >> Thanks >> >> Barry >> >> On Apr 23, 2012, at 5:50 PM, Mark F. Adams wrote: >> >>> >>> On Apr 23, 2012, at 5:59 PM, Barry Smith wrote: >>> >>>> >>>> It was changed a while ago that MatSetBlockSize() couldn't be set after the matrix was full instantiated. I guess that example did not get fixed because it is not listed in the makefile to run in the makeall. >>> >>> May I add it? >>> >>>> >>>> I have fixed the example to call MatSetBlockSize() at a safe point. >>>> >>>> Barry >>>> >>>> On Apr 23, 2012, at 4:29 PM, Mark F. Adams wrote: >>>> >>>>> ex55.c in ksp is failing with: >>>>> >>>>> [0]PETSC ERROR: --------------------- Error Message ------------------------------------ >>>>> [0]PETSC ERROR: Arguments are incompatible! >>>>> [0]PETSC ERROR: Cannot change block size 1 to 2! >>>>> [0]PETSC ERROR: ------------------------------------------------------------------------ >>>>> >>>>> on this line 57 of ex55.c: >>>>> >>>>> ierr = MatSetBlockSize(Amat,2); CHKERRQ(ierr); >>>>> >>>>> Any idea what happened here? >>>>> >>>>> Mark >>>> >>>> >>> >> >> >

1 0

Re: [petsc-dev] MetSetBlockSize issue
by Mark F. Adams 24 Apr '12

24 Apr '12

Now I'm getting this error with code like this: ierr = MatGetSubMatrix( Cmat, new_eq_indices, new_eq_indices, MAT_INITIAL_MATRIX, &mat ); CHKERRQ(ierr); ierr = MatSetBlockSize( mat, cbs ); CHKERRQ(ierr); and like this: ierr = MatPtAP( Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, &Cmat ); CHKERRQ(ierr); ierr = MatSetBlockSize( Cmat, cbs ); CHKERRQ(ierr); Mark On Apr 23, 2012, at 9:09 PM, Barry Smith wrote: > > Yes, look, for example how ex32 is handled at the bottom of the makefile. > > Thanks > > Barry > > On Apr 23, 2012, at 5:50 PM, Mark F. Adams wrote: > >> >> On Apr 23, 2012, at 5:59 PM, Barry Smith wrote: >> >>> >>> It was changed a while ago that MatSetBlockSize() couldn't be set after the matrix was full instantiated. I guess that example did not get fixed because it is not listed in the makefile to run in the makeall. >> >> May I add it? >> >>> >>> I have fixed the example to call MatSetBlockSize() at a safe point. >>> >>> Barry >>> >>> On Apr 23, 2012, at 4:29 PM, Mark F. Adams wrote: >>> >>>> ex55.c in ksp is failing with: >>>> >>>> [0]PETSC ERROR: --------------------- Error Message ------------------------------------ >>>> [0]PETSC ERROR: Arguments are incompatible! >>>> [0]PETSC ERROR: Cannot change block size 1 to 2! >>>> [0]PETSC ERROR: ------------------------------------------------------------------------ >>>> >>>> on this line 57 of ex55.c: >>>> >>>> ierr = MatSetBlockSize(Amat,2); CHKERRQ(ierr); >>>> >>>> Any idea what happened here? >>>> >>>> Mark >>> >>> >> > >

1 0

Performance of VecMDot_SeqCUSP
by Jose E. Roman 24 Apr '12

24 Apr '12

It seems that VecMDot_SeqCUSP has rather poor performance. This has a lot of impact in SLEPc because it is the main kernel used in the orthogonalization of vectors. Is this due to the version of Thrust? I am using CUDA Toolkit 4.0. I tried a naive replacement that copies the contents of the vectors into a matrix and calls CUBLAS dgemv. The improvement is significant, despite the data movement overhead. In some tests I see a reduction of time (VecReduceArith) from 24.5 seconds to 9.6 seconds (with up to 200 vectors of length 10000) on a Fermi. I can send the code for you to try. Jose

1 0

Re: [petsc-dev] MatAXPY with DIFFERENT_NONZERO_PATTERN
by Jed Brown 24 Apr '12

24 Apr '12

MPI has some asinine semantics with MPI_Ops being used differently with collectives and with one-sided. I pushed a fix yesterday. On Mon, Apr 23, 2012 at 08:46, Alexander Grayver <agrayver(a)gfz-potsdam.de>wrote: > On 23.04.2012 15:12, Jed Brown wrote: > >> My fault, now fixed: >> >> http://petsc.cs.iit.edu/petsc/**petsc-dev/rev/8a997a5a65d3<http://petsc.cs.iit.edu/petsc/petsc-dev/rev/8a997a5a65d3> >> > > Jed, > > MatAXPY works fine, although during tests there was crash: > > Possible error running C/C++ src/snes/examples/tutorials/**ex19 with 2 > MPI processes > See http://www.mcs.anl.gov/petsc/**documentation/faq.html<http://www.mcs.anl.gov/petsc/documentation/faq.html> > lid velocity = 0.0016, prandtl # = 1, grashof # = 1 > [0]PETSC ERROR: PetscSplitReductionApply() line 328 in > /home/lib/petsc-dev1/src/vec/**vec/utils/comb.c > [0]PETSC ERROR: PetscSplitReductionEnd() line 247 in > /home/lib/petsc-dev1/src/vec/**vec/utils/comb.c > [0]PETSC ERROR: VecNormEnd() line 672 in /home/lib/petsc-dev1/src/vec/** > vec/utils/comb.c > [0]PETSC ERROR: SNESSolve_LS() line 172 in /home/lib/petsc-dev1/src/snes/* > *impls/ls/ls.c > [0]PETSC ERROR: SNESSolve() line 3458 in /home/lib/petsc-dev1/src/snes/** > interface/snes.c > [0]PETSC ERROR: main() line 151 in src/snes/examples/tutorials/**ex19.c > [1]PETSC ERROR: PetscSplitReductionApply() line 328 in > /home/lib/petsc-dev1/src/vec/**vec/utils/comb.c > [1]PETSC ERROR: PetscSplitReductionEnd() line 247 in > /home/lib/petsc-dev1/src/vec/**vec/utils/comb.c > [1]PETSC ERROR: VecNormEnd() line 672 in /home/lib/petsc-dev1/src/vec/** > vec/utils/comb.c > [1]PETSC ERROR: SNESSolve_LS() line 172 in /home/lib/petsc-dev1/src/snes/* > *impls/ls/ls.c > [1]PETSC ERROR: SNESSolve() line 3458 in /home/lib/petsc-dev1/src/snes/** > interface/snes.c > [1]PETSC ERROR: main() line 151 in src/snes/examples/tutorials/**ex19.c > > -- > Regards, > Alexander > >

1 0

Re: [petsc-dev] Detecting MPI_C_DOUBLE_COMPLEX when clanguage!='C'
by Jed Brown 24 Apr '12

24 Apr '12

On Tue, Apr 24, 2012 at 09:02, Jack Poulson <jack.poulson(a)gmail.com> wrote: > On Tue, Apr 24, 2012 at 7:19 AM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > >> On Tue, Apr 24, 2012 at 07:15, Matthew Knepley <knepley(a)gmail.com> wrote: >> >>> Of course don't fix it. Who do these guys work for, Sandia? >> >> >> They could fix it as an extension in MPICH2, but it would have to be >> added to the standard for us to rely on it. Considering that we still can't >> use anything past MPI-1.1, we won't be able to remove the hacks for a long >> time. This is not to say that it isn't worth fixing, just that your >> daughter will be in college by the time the fix can be assumed by callers. >> > > Why not just cast the std::complex<double>* buffer into a a double* buffer > of twice the length and use MPI_SUM for the MPI_DOUBLE datatype? I have > always found MPI's support for complex datatypes to be hopelessly buggy, so > I play these tricks whenever possible. Is this what PETSc currently does > for complex MPI_SUM? > PETSc configure checks for MPI_C_COMPLEX_DOUBLE when using C99 complex (but not with C++ because, although the type exists in the header, it errors at run-time on Windows). For C++ complex, we do not try to use an MPI built-in (because the only way to get it is to use the deprecated C++ bindings), instead we create our own MPI_Datatype to hold complex values and our own MPI_Op that is capable of operating on complex values. This works fine for collectives, but cannot be used for MPI_Accumulate. Therefore, anywhere that we call MPI_Accumulate, we "translate" our MPI_Op to the appropriate predefined MPI_Op. After jumping through these hoops, at least the user can just always use our MPI_Ops regardless of whether the interface is implemented using collectives or one-sided.

1 0

Re: [petsc-dev] Detecting MPI_C_DOUBLE_COMPLEX when clanguage!='C'
by Jack Poulson 24 Apr '12

24 Apr '12

On Tue, Apr 24, 2012 at 7:19 AM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > On Tue, Apr 24, 2012 at 07:15, Matthew Knepley <knepley(a)gmail.com> wrote: > >> Of course don't fix it. Who do these guys work for, Sandia? > > > They could fix it as an extension in MPICH2, but it would have to be added > to the standard for us to rely on it. Considering that we still can't use > anything past MPI-1.1, we won't be able to remove the hacks for a long > time. This is not to say that it isn't worth fixing, just that your > daughter will be in college by the time the fix can be assumed by callers. > Why not just cast the std::complex<double>* buffer into a a double* buffer of twice the length and use MPI_SUM for the MPI_DOUBLE datatype? I have always found MPI's support for complex datatypes to be hopelessly buggy, so I play these tricks whenever possible. Is this what PETSc currently does for complex MPI_SUM? Jack

1 0

Re: [petsc-dev] Detecting MPI_C_DOUBLE_COMPLEX when clanguage!='C'
by Barry Smith 24 Apr '12

24 Apr '12

What purpose does it serve for you to watch a particular ticket and then see that it is now "resolved" even though the problem still exists for you. That only makes you cynical about the package and its developers. On Apr 24, 2012, at 8:33 AM, Jed Brown wrote: > The alternative is that > > 1. Doing nothing means it gets lost without resolution. > > 2. The only way to be notified of action is to subscribe to the mailing list (and not all fixes are mentioned there) which is high volume if all you care about is a particular bug. > > > On Apr 24, 2012 8:05 AM, "Barry Smith" <bsmith(a)mcs.anl.gov> wrote: > >> > >> > >> From that page "status changed from new to closed" > >> > >> This is why I absolutely hate the idea of using "ticket" systems for open source projects. The goal becomes "close the tickets at all cost using any lame ass excuse you can" instead of actually thinking about the problem and coming up with a satisfactory solution. > >> > >> BTW: How is OpenMPI on complex for Windows and C++ and stuff? > >> > >> Barry > >> > >> On Apr 24, 2012, at 7:24 AM, Matthew Knepley wrote: > >> > >> > On Tue, Apr 24, 2012 at 8:19 AM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > >> > On Tue, Apr 24, 2012 at 07:15, Matthew Knepley <knepley(a)gmail.com> wrote: > >> > Of course don't fix it. Who do these guys work for, Sandia? > >> > > >> > They could fix it as an extension in MPICH2, but it would have to be added to the standard for us to rely on it. Considering that we still can't use anything past MPI-1.1, we won't be able to remove the hacks for a long time. This is not to say that it isn't worth fixing, just that your daughter will be in college by the time the fix can be assumed by callers. > >> > > >> > MPICH used to lead by just doing good, useful things. The embarrassed standards committee > >> > would then adopt them. Now this is reversed? > >> > > >> > Matt > >> > > >> > -- > >> > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > >> > -- Norbert Wiener > >>

1 0