Re: [petsc-dev] MatHermitianTransposeGetMat

2018-10-19 Thread Pierre Jolivet
Hi Dave,

> On 19 Oct 2018, at 11:10 PM, Dave May  wrote:
> 
> Hi Pierre,
> 
> On Fri, 19 Oct 2018 at 21:28, Pierre Jolivet  > wrote:
> Hello,
> I’m trying to call MatTransposeGetMat on a Mat created with 
> MatCreateHermitianTranspose, but this fails.
> Why is there no MatHermitianTransposeGetMat?
> 
> Looks like an oversight to me.

That is what I was afraid of. I’ve proposed a pull request as you suggested.

Thanks,
Pierre

> Is there a work around?
> 
> Doesn't seem to be a hack workaround possible as the object holding the 
> reference to the original matrix is declared directly in the C file. See 
> typedef struct Mat_HT
> 
> https://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/transpose/htransm.c.html
>  
> 
> 
> Support for XXXGetMat() needs to go into master. Essentially one needs to 
> just copy the code from 
> 
> https://www.mcs.anl.gov/pets/petsc-current/src/mat/impls/transpose/transm.c.html
>  
> 
> associated with MatTransposeGetMat and change the type from Mat_Transpose to 
> Mat_HT.
> 
> Want to make a PR?
> 
> Cheers,
>   Dave
> 
> 
> 
> 
> 
> 
> 
> Thanks,
> Pierre



[petsc-dev] MatNest and MatCreateHermitianTranspose

2018-10-19 Thread Pierre Jolivet
Hello,
I’m using a MatNest with some blocks made of matrices created with 
MatCreateHermitianTranspose. I end up with the following trace when doing a 
simple MatMult.
Should this be fixed in matnest.c, matrix.c or htransm.c?

Thanks,
Pierre

(lldb) down
frame #11: 0x00011ad06830 
libpetsc.3.010.dylib`MatMult_Nest(A=0x7fcc2eb50260, x=0x7fcc2f019460, 
y=0x7fcc2f01ca60) at matnest.c:52
   49   for (j=0; jm[i][j]) continue;
   51 /* y[i] <- y[i] + A[i][j] * x[j] */
-> 52 ierr = MatMultAdd(bA->m[i][j],bx[j],by[i],by[i]);CHKERRQ(ierr);
   53   }
   54 }
   55 for (i=0; iisglobal.row[i],[i]);CHKERRQ(ierr);}
(lldb) 
frame #10: 0x00011a9d6399 
libpetsc.3.010.dylib`MatMultAdd(mat=0x7fcc2ebfde60, v1=0x7fcc2f029c60, 
v2=0x7fcc2f024660, v3=0x7fcc2f024660) at matrix.c:2517
   2514   if (!mat->ops->multadd) 
SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No MatMultAdd() for 
matrix type '%s'",((PetscObject)mat)->type_name);
   2515   ierr = PetscLogEventBegin(MAT_MultAdd,mat,v1,v2,v3);CHKERRQ(ierr);
   2516   ierr = VecLockPush(v1);CHKERRQ(ierr);
-> 2517   ierr = (*mat->ops->multadd)(mat,v1,v2,v3);CHKERRQ(ierr);
   2518   ierr = VecLockPop(v1);CHKERRQ(ierr);
   2519   ierr = PetscLogEventEnd(MAT_MultAdd,mat,v1,v2,v3);CHKERRQ(ierr);
   2520   ierr = PetscObjectStateIncrease((PetscObject)v3);CHKERRQ(ierr);
(lldb) 
frame #9: 0x00011ad136d0 
libpetsc.3.010.dylib`MatMultAdd_HT(N=0x7fcc2ebfde60, v1=0x7fcc2f029c60, 
v2=0x7fcc2f024660, v3=0x7fcc2f024660) at htransm.c:24
   21 PetscErrorCode ierr;
   22   
   23 PetscFunctionBegin;
-> 24 ierr = MatMultHermitianTransposeAdd(Na->A,v1,v2,v3);CHKERRQ(ierr);
   25 PetscFunctionReturn(0);
   26   }
   27   
(lldb) 
frame #8: 0x00011a9d88e5 
libpetsc.3.010.dylib`MatMultHermitianTransposeAdd(mat=0x7fcc2e977260, 
v1=0x7fcc2f029c60, v2=0x7fcc2f024660, v3=0x7fcc2f024660) at 
matrix.c:2629
   2626 ierr = MatMultTranspose(mat,w,z);CHKERRQ(ierr);
   2627 ierr = VecDestroy();CHKERRQ(ierr);
   2628 ierr = VecConjugate(z);CHKERRQ(ierr);
-> 2629 ierr = VecWAXPY(v3,1.0,v2,z);CHKERRQ(ierr);
   2630 ierr = VecDestroy();CHKERRQ(ierr);
   2631   }
   2632   ierr = VecLockPop(v1);CHKERRQ(ierr);
(lldb) 
frame #7: 0x00011a8f33b5 
libpetsc.3.010.dylib`VecWAXPY(w=0x7fcc2f024660, alpha=1, 
x=0x7fcc2f024660, y=0x7fcc2f039260) at rvector.c:795
   792VecCheckSameSize(x,3,y,4);
   793VecCheckSameSize(x,3,w,1);
   794if (w == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w 
cannot be same as input vector y, suggest VecAXPY()");
-> 795if (w == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w 
cannot be same as input vector x, suggest VecAYPX()");
   796PetscValidLogicalCollectiveScalar(y,alpha,2);
   797  
   798ierr = PetscLogEventBegin(VEC_WAXPY,x,y,w,0);CHKERRQ(ierr);



Re: [petsc-dev] MatHermitianTransposeGetMat

2018-10-19 Thread Dave May
Hi Pierre,

On Fri, 19 Oct 2018 at 21:28, Pierre Jolivet 
wrote:

> Hello,
> I’m trying to call MatTransposeGetMat on a Mat created with
> MatCreateHermitianTranspose, but this fails.
> Why is there no MatHermitianTransposeGetMat?


Looks like an oversight to me.


Is there a work around?


Doesn't seem to be a hack workaround possible as the object holding the
reference to the original matrix is declared directly in the C file. See
typedef struct Mat_HT

https://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/transpose/htransm.c.html

Support for XXXGetMat() needs to go into master. Essentially one needs to
just copy the code from

https://www.mcs.anl.gov/pets/petsc-current/src/mat/impls/transpose/transm.c.html

associated with MatTransposeGetMat and change the type from Mat_Transpose
to Mat_HT.

Want to make a PR?

Cheers,
  Dave






>
> Thanks,
> Pierre


Re: [petsc-dev] Fortran interface problem in v.3.10.2

2018-10-19 Thread Smith, Barry F.

  Adrian,

  I goggled a bit but couldn't understand anything I found. My guess is that

type(*)

doesn't work for certain derived types (why not?) perhaps those that contain 
procedures. If I remove the procedure from the
routine it all compiles, thus producing some evidence my theory is correct.

  So we have a problem. The type checking one user wants is causing another 
users code to not build.

  Here is a short term fix you can do.  After you run ./configure edit 
$PETSC_ARCH/include/petscconf.h locate

#ifndef PETSC_HAVE_FORTRAN_TYPE_STAR
#define PETSC_HAVE_FORTRAN_TYPE_STAR 1
#endif

then run make. This will turn off the type checking for all routines that have 
type(*) arguments which includes SNESSetConvergenceTest

  I don't have a good long term solution at the moment. Maybe a Fortran expert 
has some idea.

   Barry







> On Oct 18, 2018, at 10:14 PM, Adrian Croucher  
> wrote:
>
> hi Barry,
>
> On 18/10/18 11:34 AM, Smith, Barry F. wrote:
>>Sorry about this problem. I think the change was only introduced in 
>> master and should not affect 3.10.x Please confirm that master is where the 
>> failed compile is?
>
> Yes, it's on master. I have tracked down the commit at which it started to 
> fail: 6f222c9d1e, "type checking for Fortran" (Fri Sep 28).
>
>>Please send us the calling sequence of your routine that won't compile 
>> (cut and paste).
>
> I've attached a minimal example program which fails with the following error:
>
>call SNESSetConvergenceTest(snes, convergence, context, &
>  1
> Error: Actual argument at (1) to assumed-type dummy is of derived type with 
> type-bound or FINAL procedures
>
> Cheers, Adrian
>
>>> On Oct 17, 2018, at 5:21 PM, Adrian Croucher  
>>> wrote:
>>>
>>> hi
>>>
>>> A colleague has just reported that my code no longer builds with PETSc 
>>> 3.10.2, though it builds OK with 3.10.1.
>>>
>>> The problem appears to be the Fortran interface to 
>>> SNESSetConvergenceTest(), which was changed at commit f9a1a4d.
>>>
>>> It now complains about the context argument we are passing in to this 
>>> function ('cctx' in the interface, which is declared there as type(*)). The 
>>> error is:
>>>
>>> "Error: Actual argument at (1) to assumed-type dummy is of derived type 
>>> with type-bound or FINAL procedures"
>>>
>>> This is true, the argument being passed in is of derived type with 
>>> type-bound procedures. Previously this didn't bother it, but it looks like 
>>> it does now.
>>>
>>> Is its complaint legitimate? or perhaps a compiler bug? (this is using gcc 
>>> 6.3.0)
>>>
>>> - Adrian
>>>
> --
> Dr Adrian Croucher
> Senior Research Fellow
> Department of Engineering Science
> University of Auckland, New Zealand
> email: a.crouc...@auckland.ac.nz
> tel: +64 (0)9 923 4611
>
> 



test_snes.F90
Description: test_snes.F90


[petsc-dev] MatHermitianTransposeGetMat

2018-10-19 Thread Pierre Jolivet
Hello,
I’m trying to call MatTransposeGetMat on a Mat created with 
MatCreateHermitianTranspose, but this fails.
Why is there no MatHermitianTransposeGetMat? Is there a work around?

Thanks,
Pierre