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

Re: [petsc-dev] MatMultHermitianTranspose
by Jose E. Roman 18 Jul '12

18 Jul '12

El 18/07/2012, a las 17:45, Jed Brown escribió: > On Wed, Jul 18, 2012 at 10:22 AM, Jose E. Roman <jroman(a)dsic.upv.es> wrote: > El 18/07/2012, a las 15:05, Jed Brown escribió: > > > MatMultHermitianTranspose could have a default implementation in terms of the others. Then it would work in all cases, right? > > > > Do you mean VecConjugate(x) plus MatMultTranspose(A,x,y) plus VecConjugate(y) ?? > I guess this would work. > > In the PETSc implementation, right? Yes, I think this belongs to PETSc. > > Also, I notice this in your current code. It is always safe to call MatConjugate(), the function will just do nothing if not using complex. I see. Thanks. > > src/svd/interface/svdsetup.c- ierr = MatTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->AT);CHKERRQ(ierr); > src/svd/interface/svdsetup.c-#if defined(PETSC_USE_COMPLEX) > src/svd/interface/svdsetup.c: ierr = MatConjugate(svd->AT);CHKERRQ(ierr); > src/svd/interface/svdsetup.c-#endif > src/svd/interface/svdsetup.c- } else { > src/svd/interface/svdsetup.c- ierr = MatTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->A);CHKERRQ(ierr); > src/svd/interface/svdsetup.c-#if defined(PETSC_USE_COMPLEX) > src/svd/interface/svdsetup.c: ierr = MatConjugate(svd->A);CHKERRQ(ierr); > src/svd/interface/svdsetup.c-#endif >

1 0

Re: [petsc-dev] MatMultHermitianTranspose
by Jed Brown 18 Jul '12

18 Jul '12

On Wed, Jul 18, 2012 at 10:22 AM, Jose E. Roman <jroman(a)dsic.upv.es> wrote: > El 18/07/2012, a las 15:05, Jed Brown escribió: > > > MatMultHermitianTranspose could have a default implementation in terms > of the others. Then it would work in all cases, right? > > > > Do you mean VecConjugate(x) plus MatMultTranspose(A,x,y) plus > VecConjugate(y) ?? > I guess this would work. > In the PETSc implementation, right? Also, I notice this in your current code. It is always safe to call MatConjugate(), the function will just do nothing if not using complex. src/svd/interface/svdsetup.c- ierr = MatTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->AT);CHKERRQ(ierr); src/svd/interface/svdsetup.c-#if defined(PETSC_USE_COMPLEX) src/svd/interface/svdsetup.c: ierr = MatConjugate(svd->AT);CHKERRQ(ierr); src/svd/interface/svdsetup.c-#endif src/svd/interface/svdsetup.c- } else { src/svd/interface/svdsetup.c- ierr = MatTranspose(svd->OP,MAT_INITIAL_MATRIX,&svd->A);CHKERRQ(ierr); src/svd/interface/svdsetup.c-#if defined(PETSC_USE_COMPLEX) src/svd/interface/svdsetup.c: ierr = MatConjugate(svd->A);CHKERRQ(ierr); src/svd/interface/svdsetup.c-#endif

1 0

Re: [petsc-dev] MatMultHermitianTranspose
by Jose E. Roman 18 Jul '12

18 Jul '12

El 18/07/2012, a las 15:05, Jed Brown escribió: > MatMultHermitianTranspose could have a default implementation in terms of the others. Then it would work in all cases, right? > Do you mean VecConjugate(x) plus MatMultTranspose(A,x,y) plus VecConjugate(y) ?? I guess this would work. Jose

1 0

Re: [petsc-dev] MatMultHermitianTranspose
by Jed Brown 18 Jul '12

18 Jul '12

MatMultHermitianTranspose could have a default implementation in terms of the others. Then it would work in all cases, right? On Wed, Jul 18, 2012 at 7:59 AM, Jose E. Roman <jroman(a)dsic.upv.es> wrote: > Now the SVD in SLEPc correctly applies complex conjugation in complex > scalar builds. We have two options: > - Implicit: uses MatMultHermitianTranspose > - Explicit: uses MatTranspose plus MatConjugate (the default) > > The problem is that MatMultHermitianTranspose is defined only for SEQBAIJ > (not for parallel BAIJ nor for AIJ sequential or parallel). So the implicit > option is quite limited now for complex. > > Conversely, the explicit option works in all cases except for BAIJ > matrices since MatConjugate is not defined for BAIJ. > > How difficult is to implement these operations? > > Jose > >

1 0

MatMultHermitianTranspose
by Jose E. Roman 18 Jul '12

18 Jul '12

Now the SVD in SLEPc correctly applies complex conjugation in complex scalar builds. We have two options: - Implicit: uses MatMultHermitianTranspose - Explicit: uses MatTranspose plus MatConjugate (the default) The problem is that MatMultHermitianTranspose is defined only for SEQBAIJ (not for parallel BAIJ nor for AIJ sequential or parallel). So the implicit option is quite limited now for complex. Conversely, the explicit option works in all cases except for BAIJ matrices since MatConjugate is not defined for BAIJ. How difficult is to implement these operations? Jose

1 0

Re: [petsc-dev] where is type(Mat) defined?
by Blaise Bourdin 17 Jul '12

17 Jul '12

On Jul 17, 2012, at 10:51 PM, Barry Smith wrote: > > If you are using Fortran modules (and being a FORTRAN lover you are :-). Then you have to rebuild the modules whenever you change these files. Do make cmake after you edit the file You bet I am... make cmake is what I was missing. Is it possible to apply the attached patch to petsc-3.3 and petsc-dev? It adds the data type, and the fortran binding for MatNullSpaceView. Blaise > > Barry > > On Jul 17, 2012, at 10:41 PM, Blaise Bourdin wrote: > >> Hi, >> >> I am still trying to figure out how to add support for MatNullSpace with fortran data types. As I said in a previous email, adding the statement >> type MatNullSpace >> PetscFortranAddr:: v >> end type MatNullSpace >> in include/finclude/ftn-custom/petscmatdef.h90 >> >> is not enough. >> >> Furthermore, if I remove the corresponding statement >> type Mat >> PetscFortranAddr:: v >> end type Mat >> for Mat, I can still compile code with a statement like >> Type(Mat) :: M >> >> It looks like the fortran data type for Mat is defined somewhere else, but I can't figure out where. Can anybody shine some light for me? >> >> Blaise >> -- >> Department of Mathematics and Center for Computation & Technology >> Louisiana State University, Baton Rouge, LA 70803, USA >> Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin >> >> >> >> >> >> >> > -- Department of Mathematics and Center for Computation & Technology Louisiana State University, Baton Rouge, LA 70803, USA Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

1 0

Re: [petsc-dev] where is type(Mat) defined?
by Barry Smith 17 Jul '12

17 Jul '12

If you are using Fortran modules (and being a FORTRAN lover you are :-). Then you have to rebuild the modules whenever you change these files. Do make cmake after you edit the file Barry On Jul 17, 2012, at 10:41 PM, Blaise Bourdin wrote: > Hi, > > I am still trying to figure out how to add support for MatNullSpace with fortran data types. As I said in a previous email, adding the statement > type MatNullSpace > PetscFortranAddr:: v > end type MatNullSpace > in include/finclude/ftn-custom/petscmatdef.h90 > > is not enough. > > Furthermore, if I remove the corresponding statement > type Mat > PetscFortranAddr:: v > end type Mat > for Mat, I can still compile code with a statement like > Type(Mat) :: M > > It looks like the fortran data type for Mat is defined somewhere else, but I can't figure out where. Can anybody shine some light for me? > > Blaise > -- > Department of Mathematics and Center for Computation & Technology > Louisiana State University, Baton Rouge, LA 70803, USA > Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin > > > > > > >

1 0

Re: [petsc-dev] threadcomm memory leak
by Barry Smith 17 Jul '12

17 Jul '12

On Jul 17, 2012, at 10:14 PM, Jed Brown wrote: > On Tue, Jul 17, 2012 at 10:10 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > But Shri indicated there is a possible problem, what is that problem? Or is there none? > > The "problem" is that there is special code that removes the threadcomm from SELF and WORLD. If we assume that they are different, but in fact, SELF and WORLD are the same object, it looks like a reference counting inconsistency. Thus special code to deal with MPIUNI where they are the same. I don't understand why it's such a problem to make MPIUNI non-special in having SELF and WORLD be different values. Cause I'm a crotchety old bastard.

1 0

where is type(Mat) defined?
by Blaise Bourdin 17 Jul '12

17 Jul '12

Hi, I am still trying to figure out how to add support for MatNullSpace with fortran data types. As I said in a previous email, adding the statement type MatNullSpace PetscFortranAddr:: v end type MatNullSpace in include/finclude/ftn-custom/petscmatdef.h90 is not enough. Furthermore, if I remove the corresponding statement type Mat PetscFortranAddr:: v end type Mat for Mat, I can still compile code with a statement like Type(Mat) :: M It looks like the fortran data type for Mat is defined somewhere else, but I can't figure out where. Can anybody shine some light for me? Blaise -- Department of Mathematics and Center for Computation & Technology Louisiana State University, Baton Rouge, LA 70803, USA Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

1 0

Re: [petsc-dev] threadcomm memory leak
by Shri 17 Jul '12

17 Jul '12

On Jul 17, 2012, at 10:10 PM, Barry Smith wrote: > > On Jul 17, 2012, at 10:07 PM, Jed Brown wrote: > >> On Tue, Jul 17, 2012 at 10:06 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: >> Ok, then the destructor called by the MPI_Comm_del() would decrease the reference count by one and if it gets to zero then does the delete as usual. >> >> This is exactly what is done, but PETSC_COMM_WORLD and PETSC_COMM_SELF outlive PetscFinalize, so we remove one more reference count for each of those so that it actually gets collected. > > But Shri indicated there is a possible problem, what is that problem? Or is there none? There is no problem. I was just trying to put across a generalized way of handling MPI attributes in the PETSc code so that specific implementations, e.g threadcomm and others, can use a common code rather than having specific implementations. Shri > > Barry > >

1 0