Re: [petsc-users] With-batch (new) flags

[Smith, Barry F. via petsc-users]

  Yes, this is totally my fault. By removing the help message it made configure 
treat the argument as a string hence '0' was true and you got the error 
message. For fblaslapack one should use -known-64-bit-blas-indices=0 just as 
you did, I have pushed a fix to master

  What kind of system is sunfire09.pppl.gov ? Surely a system that has a batch 
system provides its own good BLAS/LAPACK. You should use the ones on the 
machine, not fblaslapack. Using fblaslapack in this situation is like going to 
a fancy sit-down dinner but bringing your dessert from McDonalds.

  It may be possible to remove many (but not all) of the cases where 
-known-64-bit-blas-indices is needed (for example when MKL, fblaslapack, 
f2blaslapack or --download-openblas is used we we know if the library is 64 bit 
indices and should set that without a need for a test or command line option. 
I'll look at it.

  Barry


> On May 20, 2019, at 3:50 PM, Balay, Satish via petsc-users 
>  wrote:
> 
> I'm not yet sure what the correct fix is - but the following change should 
> get this going..
> 
> diff --git a/config/BuildSystem/config/packages/BlasLapack.py 
> b/config/BuildSystem/config/packages/BlasLapack.py
> index e0310da4b0..7355f1a369 100644
> --- a/config/BuildSystem/config/packages/BlasLapack.py
> +++ b/config/BuildSystem/config/packages/BlasLapack.py
> @@ -42,7 +42,7 @@ class Configure(config.package.Package):
> help.addArgument('BLAS/LAPACK', '-with-lapack-lib= [/Users//liblapack.a,...]>',nargs.ArgLibrary(None, None, 'Indicate the 
> library(s) containing LAPACK'))
> help.addArgument('BLAS/LAPACK', 
> '-with-blaslapack-suffix=',nargs.ArgLibrary(None, None, 'Indicate a 
> suffix for BLAS/LAPACK subroutine names.'))
> help.addArgument('BLAS/LAPACK', '-with-64-bit-blas-indices', 
> nargs.ArgBool(None, 0, 'Try to use 64 bit integers for BLAS/LAPACK; will 
> error if not available'))
> -#help.addArgument('BLAS/LAPACK', '-known-64-bit-blas-indices=', 
> nargs.ArgBool(None, 0, 'Indicate if using 64 bit integer BLAS'))
> +help.addArgument('BLAS/LAPACK', '-known-64-bit-blas-indices=', 
> nargs.ArgBool(None, 0, 'Indicate if using 64 bit integer BLAS'))
> return
> 
>   def getPrefix(self):
> 
> Satish
> 
> On Mon, 20 May 2019, Mark Adams via petsc-users wrote:
> 
>> On Mon, May 20, 2019 at 3:55 PM Balay, Satish  wrote:
>> 
>>> for ex:  ilp version of mkl is --known-64-bit-blas-indices=1 while lp mkl
>>> is --known-64-bit-blas-indices=0
>>> 
>>> Default blas we normally use is --known-64-bit-blas-indices=0 [they don't
>>> use 64bit indices]
>>> 
>> 
>> Humm, that is what Dylan (in the log that I sent). He is downloading blas
>> and has --known-64-bit-blas-indices=0. Should this be correct?
>> 
>> 
>>> 
>>> Satish
>>> 
>>> On Mon, 20 May 2019, Mark Adams via petsc-users wrote:
>>> 
 We are getting this failure. This a bit frustrating in that the first
>>> error
 message "Must give a default value for known-mpi-shared-libraries.." OK,
>>> I
 google it and find that =0 is suggested. That seemed to work. Then we
>>> got a
 similar error about -known-64-bit-blas-indices. It was clear from the
 documentation what to use so we tried =0 and that failed (attached). This
 is little frustrating having to use try and error for each of these
>>> 'known"
 things.
 
 Dylan is trying  --known-64-bit-blas-indices=1 now. I trust that will
>>> work,
 but I think the error are not very informative. All this known stuff is
>>> new
 to me. Perhaps put an FAQ for this and list all of the "known"s that we
 need to add in batch.
 
 Thanks,
 Mark
 
>>> 
>>> 
>> 
> 

Re: [petsc-users] Question about TSComputeRHSJacobianConstant

[Zhang, Hong via petsc-users]

Sajid,

I have also rested the simpler problem you provided. The branch 
hongzh/fix-computejacobian gives exactly the same numerical results as the 
master branch does, but runs much faster. So the solver seems to work correctly.

To rule out the possible compiler issues, you might want to try a different 
compiler or different optimization flags. Also you might want to try smaller 
stepsizes.

Hong

On May 20, 2019, at 4:47 PM, Sajid Ali 
mailto:sajidsyed2...@u.northwestern.edu>> 
wrote:

Hi Hong,

I tried running a simpler problem that solves the equation ` u_t = A*u_xx + 
A*u_yy; ` and the fix-computejacobian branch works for this on a coarse grid. 
The code for the same is here : 
https://github.com/s-sajid-ali/xwp_petsc/blob/master/2d/FD/free_space/ex_dmda.c 
and it requires no input file. It writes at all times steps and comparing the 
output at the last time step, everything looks fine.

I want to eliminate a possible source of error which could be the fact that I 
installed both versions (3.11.2 and fix-computejacobian) with intel compilers 
and O3. Could floating point errors occur due to this ? I didn't specify 
-fp-model strict but since the results I got were reasonable I never bothered 
to run a test suite.

Thank You,
Sajid Ali
Applied Physics
Northwestern University

Re: [petsc-users] DMNetwork in petsc4py

[Smith, Barry F. via petsc-users]

Justin,

That would be great. No one is working on it that I know of.

   Barry

> On May 20, 2019, at 1:48 PM, Justin Chang via petsc-users 
>  wrote:
> 
> Hi all,
> 
> Is there any current effort or plan to make the DMNetwork calls available in 
> petsc4py? I don't see anything DMNetwork related in the master branch. 
> Because we (NREL) have lots of existing network-like applications written in 
> Python and Julia, none of which use PETSc, and instead use unscalable direct 
> methods, so I think instead of rewriting everything into a C-based PETSc code 
> it would be much more feasible to integrate their application code into 
> petsc4py with DMNetwork. 
> 
> If not I was planning on writing these wrappers myself but wasn't sure if 
> someone here is already doing it.
> 
> Thanks,
> Justin

Re: [petsc-users] With-batch (new) flags

[Balay, Satish via petsc-users]

I'm not yet sure what the correct fix is - but the following change should get 
this going..

diff --git a/config/BuildSystem/config/packages/BlasLapack.py 
b/config/BuildSystem/config/packages/BlasLapack.py
index e0310da4b0..7355f1a369 100644
--- a/config/BuildSystem/config/packages/BlasLapack.py
+++ b/config/BuildSystem/config/packages/BlasLapack.py
@@ -42,7 +42,7 @@ class Configure(config.package.Package):
 help.addArgument('BLAS/LAPACK', '-with-lapack-lib=',nargs.ArgLibrary(None, None, 'Indicate the 
library(s) containing LAPACK'))
 help.addArgument('BLAS/LAPACK', 
'-with-blaslapack-suffix=',nargs.ArgLibrary(None, None, 'Indicate a 
suffix for BLAS/LAPACK subroutine names.'))
 help.addArgument('BLAS/LAPACK', '-with-64-bit-blas-indices', 
nargs.ArgBool(None, 0, 'Try to use 64 bit integers for BLAS/LAPACK; will error 
if not available'))
-#help.addArgument('BLAS/LAPACK', '-known-64-bit-blas-indices=', 
nargs.ArgBool(None, 0, 'Indicate if using 64 bit integer BLAS'))
+help.addArgument('BLAS/LAPACK', '-known-64-bit-blas-indices=', 
nargs.ArgBool(None, 0, 'Indicate if using 64 bit integer BLAS'))
 return
 
   def getPrefix(self):

Satish

On Mon, 20 May 2019, Mark Adams via petsc-users wrote:

> On Mon, May 20, 2019 at 3:55 PM Balay, Satish  wrote:
> 
> > for ex:  ilp version of mkl is --known-64-bit-blas-indices=1 while lp mkl
> > is --known-64-bit-blas-indices=0
> >
> > Default blas we normally use is --known-64-bit-blas-indices=0 [they don't
> > use 64bit indices]
> >
> 
> Humm, that is what Dylan (in the log that I sent). He is downloading blas
> and has --known-64-bit-blas-indices=0. Should this be correct?
> 
> 
> >
> > Satish
> >
> > On Mon, 20 May 2019, Mark Adams via petsc-users wrote:
> >
> > > We are getting this failure. This a bit frustrating in that the first
> > error
> > > message "Must give a default value for known-mpi-shared-libraries.." OK,
> > I
> > > google it and find that =0 is suggested. That seemed to work. Then we
> > got a
> > > similar error about -known-64-bit-blas-indices. It was clear from the
> > > documentation what to use so we tried =0 and that failed (attached). This
> > > is little frustrating having to use try and error for each of these
> > 'known"
> > > things.
> > >
> > > Dylan is trying  --known-64-bit-blas-indices=1 now. I trust that will
> > work,
> > > but I think the error are not very informative. All this known stuff is
> > new
> > > to me. Perhaps put an FAQ for this and list all of the "known"s that we
> > > need to add in batch.
> > >
> > > Thanks,
> > > Mark
> > >
> >
> >
> 

Re: [petsc-users] With-batch (new) flags

[Mark Adams via petsc-users]

On Mon, May 20, 2019 at 3:55 PM Balay, Satish  wrote:

> for ex:  ilp version of mkl is --known-64-bit-blas-indices=1 while lp mkl
> is --known-64-bit-blas-indices=0
>
> Default blas we normally use is --known-64-bit-blas-indices=0 [they don't
> use 64bit indices]
>

Humm, that is what Dylan (in the log that I sent). He is downloading blas
and has --known-64-bit-blas-indices=0. Should this be correct?


>
> Satish
>
> On Mon, 20 May 2019, Mark Adams via petsc-users wrote:
>
> > We are getting this failure. This a bit frustrating in that the first
> error
> > message "Must give a default value for known-mpi-shared-libraries.." OK,
> I
> > google it and find that =0 is suggested. That seemed to work. Then we
> got a
> > similar error about -known-64-bit-blas-indices. It was clear from the
> > documentation what to use so we tried =0 and that failed (attached). This
> > is little frustrating having to use try and error for each of these
> 'known"
> > things.
> >
> > Dylan is trying  --known-64-bit-blas-indices=1 now. I trust that will
> work,
> > but I think the error are not very informative. All this known stuff is
> new
> > to me. Perhaps put an FAQ for this and list all of the "known"s that we
> > need to add in batch.
> >
> > Thanks,
> > Mark
> >
>
>

Re: [petsc-users] With-batch (new) flags

[Balay, Satish via petsc-users]

for ex:  ilp version of mkl is --known-64-bit-blas-indices=1 while lp mkl is 
--known-64-bit-blas-indices=0

Default blas we normally use is --known-64-bit-blas-indices=0 [they don't use 
64bit indices]

Satish

On Mon, 20 May 2019, Mark Adams via petsc-users wrote:

> We are getting this failure. This a bit frustrating in that the first error
> message "Must give a default value for known-mpi-shared-libraries.." OK, I
> google it and find that =0 is suggested. That seemed to work. Then we got a
> similar error about -known-64-bit-blas-indices. It was clear from the
> documentation what to use so we tried =0 and that failed (attached). This
> is little frustrating having to use try and error for each of these 'known"
> things.
> 
> Dylan is trying  --known-64-bit-blas-indices=1 now. I trust that will work,
> but I think the error are not very informative. All this known stuff is new
> to me. Perhaps put an FAQ for this and list all of the "known"s that we
> need to add in batch.
> 
> Thanks,
> Mark
> 

Re: [petsc-users] Creating a DMNetwork from a DMPlex

[Swarnava Ghosh via petsc-users]

Hi Barry and Matt,

Maybe try building by hand in a DMNetwork using a handrawn mesh with just a
few vertices and endless and see if what you want to do makes sense
> Okay, will try to do that. Do you have any DMNetwork example which I
could follow.

I think DMNetwork is not buying you anything here. It seems to make more
sense to do it directly in Plex.
You can easily lay down a P1 element for each field so that you can
interpolate wherever you want.
> Okay, then will it be possible to do vertex partitioning with plex?
Essentially two processes can share an element but not vertex.

I would start from a clean example, such as SNES ex17. That solves
elasticity, so it has multiple fields and FEM.
The change is that you don't want to use any of the assembly functions, so
you keep the code that does data layout
and FEM discretization, but it ignore the residual/Jacobian stuff. Feel
free to ask about using the lower-level
interpolation stuff which is not as documented.
>Thanks for pointing out the reference. Could you please share the
functions for interpolation?

Sincerely,
SG


On Mon, May 20, 2019 at 4:02 AM Matthew Knepley  wrote:

> On Mon, May 20, 2019 at 3:05 AM Swarnava Ghosh via petsc-users <
> petsc-users@mcs.anl.gov> wrote:
>
>> Hi Barry,
>>
>> Thank you for your email. My planned discretization is based on the fact
>> that I need a distributed unstructured mesh, where at each vertex point I
>> perform local calculations. For these calculations, I do NOT need need to
>> assemble any global matrix. I will have fields defined at the vertices, and
>> using linear interpolation, I am planing to find the values of these fields
>> at some spatial points with are within a ball around each vertex. Once the
>> values of these fields are known within the compact support around each
>> vertex, I do local computations to calculate my unknown field. My reason
>> for having the a mesh is essentially to 1) define fields at the vertices
>> and 2) perform linear interpolation (using finite elements) at some spatial
>> points. Also the local computations around at each vertex is
>> computationally the most expensive step. In that case, having a cell
>> partitioning will result in vertices being shared among processes, which
>> will result in redundant computations.
>>
>> My idea is therefore to have DMNetwork to distribute vertices across
>> processes and use finite elements for the linear interpolation part.
>>
>
> I think DMNetwork is not buying you anything here. It seems to make more
> sense to do it directly in Plex.
> You can easily lay down a P1 element for each field so that you can
> interpolate wherever you want.
>
> I would start from a clean example, such as SNES ex17. That solves
> elasticity, so it has multiple fields and FEM.
> The change is that you don't want to use any of the assembly functions, so
> you keep the code that does data layout
> and FEM discretization, but it ignore the residual/Jacobian stuff. Feel
> free to ask about using the lower-level
> interpolation stuff which is not as documented.
>
>   Thanks,
>
>  Matt
>
>
>> Thanks,
>> SG
>>
>>
>>
>> On Sun, May 19, 2019 at 6:54 PM Smith, Barry F. 
>> wrote:
>>
>>>
>>>   I am not sure you want DMNetwork, DMNetwork has no geometry; it only
>>> has vertices and edges. Vertices are connected to other vertices through
>>> the edges. For example I can't see how one would do vertex centered finite
>>> volume methods with DMNetwork. Maybe if you said something more about your
>>> planned discretization we could figure something out.
>>>
>>> > On May 19, 2019, at 8:32 PM, Swarnava Ghosh 
>>> wrote:
>>> >
>>> > Hi Barry,
>>> >
>>> > No, the gmesh file contains a mesh and not a graph/network.
>>> > In that case, is it possible to create a DMNetwork first from the
>>> DMPlex and then distribute the DMNetwork.
>>> >
>>> > I have this case, because I want a vertex partitioning of my mesh.
>>> Domain decomposition of DMPlex gives me cell partitioning. Essentially what
>>> I want is that no two processes can share a vertex BUT that can share an
>>> edge. Similar to how a DMDA is distributed.
>>> >
>>> > Thanks,
>>> > Swarnava
>>> >
>>> > On Sun, May 19, 2019 at 4:50 PM Smith, Barry F. 
>>> wrote:
>>> >
>>> >This use case never occurred to us. Is the gmesh file containing a
>>> graph/network (as opposed to a mesh)? There seem two choices
>>> >
>>> > 1) if the gmesh file contains a graph/network one could write a gmesh
>>> reader for that case that reads directly for and constructs a DMNetwork or
>>> >
>>> > 2) write a converter for a DMPlex to DMNetwork.
>>> >
>>> >I lean toward the first
>>> >
>>> >Either way you need to understand the documentation for DMNetwork
>>> and how to build one up.
>>> >
>>> >
>>> >Barry
>>> >
>>> >
>>> > > On May 19, 2019, at 6:34 PM, Swarnava Ghosh via petsc-users <
>>> petsc-users@mcs.anl.gov> wrote:
>>> > >
>>> > > Hi Petsc users and developers,
>>> > >
>>> > > I am trying to find a way of creating a 

Re: [petsc-users] problem with generating simplicies mesh

[Stefano Zampini via petsc-users]

Matt,

The code is actually for 2d.

> On May 20, 2019, at 12:54 PM, Matthew Knepley via petsc-users 
>  wrote:
> 
> On Sun, May 19, 2019 at 9:22 AM 陳鳴諭 via petsc-users  > wrote:
> I have problem with generating simplicies mesh.
> I do as the description in DMPlexCreateBoxmesh says, but still meet error.
> 
> Stefano is right that you will need a mesh generator for a simplex mesh. 
> However, you are asking for
> a 1D mesh for which there are no generators. Since in 1D simplces and tensor 
> cells are the same, 
> just change it to tensor and it will work.
> 
>   Thanks,
> 
>  Matt
>  
> The following is the error message:
> [0]PETSC ERROR: - Error Message 
> --
> [0]PETSC ERROR: Argument out of range
> [0]PETSC ERROR: No grid generator of dimension 1 registered
> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html 
>  for trouble shooting.
> [0]PETSC ERROR: Petsc Development GIT revision: v3.11.1-723-g96d64d1  GIT 
> Date: 2019-05-15 13:23:17 +
> [0]PETSC ERROR: ./membrane on a arch-linux2-c-debug named 
> simon-System-Product-Name by simon Sun May 19 20:54:54 2019
> [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ 
> --with-fc=gfortran --download-mpich --download-fblaslapack
> [0]PETSC ERROR: #1 DMPlexGenerate() line 181 in 
> /home/simon/petsc/src/dm/impls/plex/plexgenerate.c
> [0]PETSC ERROR: #2 DMPlexCreateBoxMesh_Simplex_Internal() line 536 in 
> /home/simon/petsc/src/dm/impls/plex/plexcreate.c
> [0]PETSC ERROR: #3 DMPlexCreateBoxMesh() line 1071 in 
> /home/simon/petsc/src/dm/impls/plex/plexcreate.c
> [0]PETSC ERROR: #4 main() line 54 in /home/simon/Downloads/membrane.c
> [0]PETSC ERROR: PETSc Option Table entries:
> [0]PETSC ERROR: -dm_view
> [0]PETSC ERROR: End of Error Message ---send entire error 
> message to petsc-ma...@mcs.anl.gov--
> application called MPI_Abort(MPI_COMM_WORLD, 63) - process 0
> [unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=63
> :
> system msg for write_line failure : Bad file descriptor
> 
> 
> 
> I need some help about this, please.
> 
> Simon
> 
> 
> -- 
> 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
> 
> https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] problem with generating simplicies mesh

[Matthew Knepley via petsc-users]

On Sun, May 19, 2019 at 9:22 AM 陳鳴諭 via petsc-users 
wrote:

> I have problem with generating simplicies mesh.
> I do as the description in DMPlexCreateBoxmesh says, but still meet error.
>

Stefano is right that you will need a mesh generator for a simplex mesh.
However, you are asking for
a 1D mesh for which there are no generators. Since in 1D simplces and
tensor cells are the same,
just change it to tensor and it will work.

  Thanks,

 Matt


> The following is the error message:
> [0]PETSC ERROR: - Error Message
> --
> [0]PETSC ERROR: Argument out of range
> [0]PETSC ERROR: No grid generator of dimension 1 registered
> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
> for trouble shooting.
> [0]PETSC ERROR: Petsc Development GIT revision: v3.11.1-723-g96d64d1  GIT
> Date: 2019-05-15 13:23:17 +
> [0]PETSC ERROR: ./membrane on a arch-linux2-c-debug named
> simon-System-Product-Name by simon Sun May 19 20:54:54 2019
> [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++
> --with-fc=gfortran --download-mpich --download-fblaslapack
> [0]PETSC ERROR: #1 DMPlexGenerate() line 181 in
> /home/simon/petsc/src/dm/impls/plex/plexgenerate.c
> [0]PETSC ERROR: #2 DMPlexCreateBoxMesh_Simplex_Internal() line 536 in
> /home/simon/petsc/src/dm/impls/plex/plexcreate.c
> [0]PETSC ERROR: #3 DMPlexCreateBoxMesh() line 1071 in
> /home/simon/petsc/src/dm/impls/plex/plexcreate.c
> [0]PETSC ERROR: #4 main() line 54 in /home/simon/Downloads/membrane.c
> [0]PETSC ERROR: PETSc Option Table entries:
> [0]PETSC ERROR: -dm_view
> [0]PETSC ERROR: End of Error Message ---send entire
> error message to petsc-ma...@mcs.anl.gov--
> application called MPI_Abort(MPI_COMM_WORLD, 63) - process 0
> [unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=63
> :
> system msg for write_line failure : Bad file descriptor
>
>
>
> I need some help about this, please.
>
> Simon
>


-- 
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

https://www.cse.buffalo.edu/~knepley/ 

Re: [petsc-users] Calling LAPACK routines from PETSc

[Dave Lee via petsc-users]

Thanks Barry,

I found some helpful examples on the intel lapack site - moral of the
story: using C ordering for input matrix, but transposed output matrices
leads to a consistent solution.

Cheers, Dave.

On Mon, May 20, 2019 at 6:07 PM Smith, Barry F.  wrote:

>
>
> > On May 20, 2019, at 2:28 AM, Dave Lee  wrote:
> >
> > Thanks Jed and Barry,
> >
> > So, just to confirm,
> >
> > -- From the KSP_GMRES structure, if I call *HH(a,b), that will return
> the row a, column b entry of the Hessenberg matrix (while the back end
> array *hh_origin array is ordering using the Fortran convention)
> >
> > -- Matrices are passed into and returned from
> PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_() using Fortran indexing, and
> need to be transposed to get back to C ordering
>
>   In general, I guess depending on what you want to do with them you don'
> need to transpose them. Why would you want to? Just leave them as little
> column oriented blogs  and with them what you need directly.
>
>Just do stuff and you'll find it works out.
>
> >
> > Are both of these statements correct?
> >
> > Cheers, Dave.
> >
> > On Mon, May 20, 2019 at 4:34 PM Smith, Barry F. 
> wrote:
> >
> >The little work arrays in GMRES tend to be stored in Fortran
> ordering; there is no C style p[][] indexing into such arrays. Thus the
> arrays can safely be sent to LAPACK.  The only trick is knowing the two
> dimensions and as Jed say the "leading dimension parameter. He gave you a
> place to look
> >
> > > On May 20, 2019, at 1:24 AM, Jed Brown via petsc-users <
> petsc-users@mcs.anl.gov> wrote:
> > >
> > > Dave Lee via petsc-users  writes:
> > >
> > >> Hi Petsc,
> > >>
> > >> I'm attempting to implement a "hookstep" for the SNES trust region
> solver.
> > >> Essentially what I'm trying to do is replace the solution of the least
> > >> squares problem at the end of each GMRES solve with a modified
> solution
> > >> with a norm that is constrained to be within the size of the trust
> region.
> > >>
> > >> In order to do this I need to perform an SVD on the Hessenberg matrix,
> > >> which copying the function KSPComputeExtremeSingularValues(), I'm
> trying to
> > >> do by accessing the LAPACK function dgesvd() via the
> PetscStackCallBLAS()
> > >> machinery. One thing I'm confused about however is the ordering of
> the 2D
> > >> arrays into and out of this function, given that that C and FORTRAN
> arrays
> > >> use reverse indexing, ie: C[j+1][i+1] = F[i,j].
> > >>
> > >> Given that the Hessenberg matrix has k+1 rows and k columns, should I
> be
> > >> still be initializing this as H[row][col] and passing this into
> > >> PetscStackCallBLAS("LAPACKgesvd",LAPACKgrsvd_(...))
> > >> or should I be transposing this before passing it in?
> > >
> > > LAPACK terminology is with respect to Fortran ordering.  There is a
> > > "leading dimension" parameter so that you can operate on non-contiguous
> > > blocks.  See KSPComputeExtremeSingularValues_GMRES for an example.
> > >
> > >> Also for the left and right singular vector matrices that are
> returned by
> > >> this function, should I be transposing these before I interpret them
> as C
> > >> arrays?
> > >>
> > >> I've attached my modified version of gmres.c in case this is helpful.
> If
> > >> you grep for DRL (my initials) then you'll see my changes to the code.
> > >>
> > >> Cheers, Dave.
> > >>
> > >> /*
> > >>This file implements GMRES (a Generalized Minimal Residual) method.
> > >>Reference:  Saad and Schultz, 1986.
> > >>
> > >>
> > >>Some comments on left vs. right preconditioning, and restarts.
> > >>Left and right preconditioning.
> > >>If right preconditioning is chosen, then the problem being solved
> > >>by gmres is actually
> > >>   My =  AB^-1 y = f
> > >>so the initial residual is
> > >>  r = f - Mx
> > >>Note that B^-1 y = x or y = B x, and if x is non-zero, the initial
> > >>residual is
> > >>  r = f - A x
> > >>The final solution is then
> > >>  x = B^-1 y
> > >>
> > >>If left preconditioning is chosen, then the problem being solved is
> > >>   My = B^-1 A x = B^-1 f,
> > >>and the initial residual is
> > >>   r  = B^-1(f - Ax)
> > >>
> > >>Restarts:  Restarts are basically solves with x0 not equal to zero.
> > >>Note that we can eliminate an extra application of B^-1 between
> > >>restarts as long as we don't require that the solution at the end
> > >>of an unsuccessful gmres iteration always be the solution x.
> > >> */
> > >>
> > >> #include <../src/ksp/ksp/impls/gmres/gmresimpl.h>   /*I
> "petscksp.h"  I*/
> > >> #include  // DRL
> > >> #define GMRES_DELTA_DIRECTIONS 10
> > >> #define GMRES_DEFAULT_MAXK 30
> > >> static PetscErrorCode
> KSPGMRESUpdateHessenberg(KSP,PetscInt,PetscBool,PetscReal*);
> > >> static PetscErrorCode
> KSPGMRESBuildSoln(PetscScalar*,Vec,Vec,KSP,PetscInt);
> > >>
> > >> PetscErrorCodeKSPSetUp_GMRES(KSP ksp)
> > >> {
> > >>  PetscInt   

Re: [petsc-users] Calling LAPACK routines from PETSc

[Smith, Barry F. via petsc-users]

> On May 20, 2019, at 2:28 AM, Dave Lee  wrote:
> 
> Thanks Jed and Barry,
> 
> So, just to confirm, 
> 
> -- From the KSP_GMRES structure, if I call *HH(a,b), that will return the row 
> a, column b entry of the Hessenberg matrix (while the back end array 
> *hh_origin array is ordering using the Fortran convention)
> 
> -- Matrices are passed into and returned from 
> PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_() using Fortran indexing, and 
> need to be transposed to get back to C ordering

  In general, I guess depending on what you want to do with them you don' need 
to transpose them. Why would you want to? Just leave them as little column 
oriented blogs  and with them what you need directly.

   Just do stuff and you'll find it works out.

> 
> Are both of these statements correct?
> 
> Cheers, Dave.
> 
> On Mon, May 20, 2019 at 4:34 PM Smith, Barry F.  wrote:
> 
>The little work arrays in GMRES tend to be stored in Fortran ordering; 
> there is no C style p[][] indexing into such arrays. Thus the arrays can 
> safely be sent to LAPACK.  The only trick is knowing the two dimensions and 
> as Jed say the "leading dimension parameter. He gave you a place to look
> 
> > On May 20, 2019, at 1:24 AM, Jed Brown via petsc-users 
> >  wrote:
> > 
> > Dave Lee via petsc-users  writes:
> > 
> >> Hi Petsc,
> >> 
> >> I'm attempting to implement a "hookstep" for the SNES trust region solver.
> >> Essentially what I'm trying to do is replace the solution of the least
> >> squares problem at the end of each GMRES solve with a modified solution
> >> with a norm that is constrained to be within the size of the trust region.
> >> 
> >> In order to do this I need to perform an SVD on the Hessenberg matrix,
> >> which copying the function KSPComputeExtremeSingularValues(), I'm trying to
> >> do by accessing the LAPACK function dgesvd() via the PetscStackCallBLAS()
> >> machinery. One thing I'm confused about however is the ordering of the 2D
> >> arrays into and out of this function, given that that C and FORTRAN arrays
> >> use reverse indexing, ie: C[j+1][i+1] = F[i,j].
> >> 
> >> Given that the Hessenberg matrix has k+1 rows and k columns, should I be
> >> still be initializing this as H[row][col] and passing this into
> >> PetscStackCallBLAS("LAPACKgesvd",LAPACKgrsvd_(...))
> >> or should I be transposing this before passing it in?
> > 
> > LAPACK terminology is with respect to Fortran ordering.  There is a
> > "leading dimension" parameter so that you can operate on non-contiguous
> > blocks.  See KSPComputeExtremeSingularValues_GMRES for an example.
> > 
> >> Also for the left and right singular vector matrices that are returned by
> >> this function, should I be transposing these before I interpret them as C
> >> arrays?
> >> 
> >> I've attached my modified version of gmres.c in case this is helpful. If
> >> you grep for DRL (my initials) then you'll see my changes to the code.
> >> 
> >> Cheers, Dave.
> >> 
> >> /*
> >>This file implements GMRES (a Generalized Minimal Residual) method.
> >>Reference:  Saad and Schultz, 1986.
> >> 
> >> 
> >>Some comments on left vs. right preconditioning, and restarts.
> >>Left and right preconditioning.
> >>If right preconditioning is chosen, then the problem being solved
> >>by gmres is actually
> >>   My =  AB^-1 y = f
> >>so the initial residual is
> >>  r = f - Mx
> >>Note that B^-1 y = x or y = B x, and if x is non-zero, the initial
> >>residual is
> >>  r = f - A x
> >>The final solution is then
> >>  x = B^-1 y
> >> 
> >>If left preconditioning is chosen, then the problem being solved is
> >>   My = B^-1 A x = B^-1 f,
> >>and the initial residual is
> >>   r  = B^-1(f - Ax)
> >> 
> >>Restarts:  Restarts are basically solves with x0 not equal to zero.
> >>Note that we can eliminate an extra application of B^-1 between
> >>restarts as long as we don't require that the solution at the end
> >>of an unsuccessful gmres iteration always be the solution x.
> >> */
> >> 
> >> #include <../src/ksp/ksp/impls/gmres/gmresimpl.h>   /*I  "petscksp.h"  
> >> I*/
> >> #include  // DRL
> >> #define GMRES_DELTA_DIRECTIONS 10
> >> #define GMRES_DEFAULT_MAXK 30
> >> static PetscErrorCode 
> >> KSPGMRESUpdateHessenberg(KSP,PetscInt,PetscBool,PetscReal*);
> >> static PetscErrorCode KSPGMRESBuildSoln(PetscScalar*,Vec,Vec,KSP,PetscInt);
> >> 
> >> PetscErrorCodeKSPSetUp_GMRES(KSP ksp)
> >> {
> >>  PetscInt   hh,hes,rs,cc;
> >>  PetscErrorCode ierr;
> >>  PetscInt   max_k,k;
> >>  KSP_GMRES  *gmres = (KSP_GMRES*)ksp->data;
> >> 
> >>  PetscFunctionBegin;
> >>  max_k = gmres->max_k;  /* restart size */
> >>  hh= (max_k + 2) * (max_k + 1);
> >>  hes   = (max_k + 1) * (max_k + 1);
> >>  rs= (max_k + 2);
> >>  cc= (max_k + 1);
> >> 
> >>  ierr = 
> >> 

Re: [petsc-users] Creating a DMNetwork from a DMPlex

[Smith, Barry F. via petsc-users]

  Maybe try building by hand in a DMNetwork using a handrawn mesh with just a 
few vertices and endless and see if what you want to do makes sense


> On May 20, 2019, at 2:04 AM, Swarnava Ghosh  wrote:
> 
> Hi Barry,
> 
> Thank you for your email. My planned discretization is based on the fact that 
> I need a distributed unstructured mesh, where at each vertex point I perform 
> local calculations. For these calculations, I do NOT need need to assemble 
> any global matrix. I will have fields defined at the vertices, and using 
> linear interpolation, I am planing to find the values of these fields at some 
> spatial points with are within a ball around each vertex. Once the values of 
> these fields are known within the compact support around each vertex, I do 
> local computations to calculate my unknown field. My reason for having the a 
> mesh is essentially to 1) define fields at the vertices and 2) perform linear 
> interpolation (using finite elements) at some spatial points. Also the local 
> computations around at each vertex is computationally the most expensive 
> step. In that case, having a cell partitioning will result in vertices being 
> shared among processes, which will result in redundant computations. 
> 
> My idea is therefore to have DMNetwork to distribute vertices across 
> processes and use finite elements for the linear interpolation part.
> 
> Thanks,
> SG
> 
> 
> 
> On Sun, May 19, 2019 at 6:54 PM Smith, Barry F.  wrote:
> 
>   I am not sure you want DMNetwork, DMNetwork has no geometry; it only has 
> vertices and edges. Vertices are connected to other vertices through the 
> edges. For example I can't see how one would do vertex centered finite volume 
> methods with DMNetwork. Maybe if you said something more about your planned 
> discretization we could figure something out.
> 
> > On May 19, 2019, at 8:32 PM, Swarnava Ghosh  wrote:
> > 
> > Hi Barry,
> > 
> > No, the gmesh file contains a mesh and not a graph/network.
> > In that case, is it possible to create a DMNetwork first from the DMPlex 
> > and then distribute the DMNetwork.
> > 
> > I have this case, because I want a vertex partitioning of my mesh. Domain 
> > decomposition of DMPlex gives me cell partitioning. Essentially what I want 
> > is that no two processes can share a vertex BUT that can share an edge. 
> > Similar to how a DMDA is distributed.
> > 
> > Thanks,
> > Swarnava
> > 
> > On Sun, May 19, 2019 at 4:50 PM Smith, Barry F.  wrote:
> > 
> >This use case never occurred to us. Is the gmesh file containing a 
> > graph/network (as opposed to a mesh)? There seem two choices
> > 
> > 1) if the gmesh file contains a graph/network one could write a gmesh 
> > reader for that case that reads directly for and constructs a DMNetwork or
> > 
> > 2) write a converter for a DMPlex to DMNetwork. 
> > 
> >I lean toward the first
> > 
> >Either way you need to understand the documentation for DMNetwork and 
> > how to build one up.
> > 
> > 
> >Barry
> > 
> > 
> > > On May 19, 2019, at 6:34 PM, Swarnava Ghosh via petsc-users 
> > >  wrote:
> > > 
> > > Hi Petsc users and developers,
> > > 
> > > I am trying to find a way of creating a DMNetwork from a DMPlex. I have 
> > > read the DMPlex from a gmesh file and have it distributed.
> > > 
> > > Thanks,
> > > SG
> > 
> 

Re: [petsc-users] Calling LAPACK routines from PETSc

[Smith, Barry F. via petsc-users]

   The little work arrays in GMRES tend to be stored in Fortran ordering; there 
is no C style p[][] indexing into such arrays. Thus the arrays can safely be 
sent to LAPACK.  The only trick is knowing the two dimensions and as Jed say 
the "leading dimension parameter. He gave you a place to look

> On May 20, 2019, at 1:24 AM, Jed Brown via petsc-users 
>  wrote:
> 
> Dave Lee via petsc-users  writes:
> 
>> Hi Petsc,
>> 
>> I'm attempting to implement a "hookstep" for the SNES trust region solver.
>> Essentially what I'm trying to do is replace the solution of the least
>> squares problem at the end of each GMRES solve with a modified solution
>> with a norm that is constrained to be within the size of the trust region.
>> 
>> In order to do this I need to perform an SVD on the Hessenberg matrix,
>> which copying the function KSPComputeExtremeSingularValues(), I'm trying to
>> do by accessing the LAPACK function dgesvd() via the PetscStackCallBLAS()
>> machinery. One thing I'm confused about however is the ordering of the 2D
>> arrays into and out of this function, given that that C and FORTRAN arrays
>> use reverse indexing, ie: C[j+1][i+1] = F[i,j].
>> 
>> Given that the Hessenberg matrix has k+1 rows and k columns, should I be
>> still be initializing this as H[row][col] and passing this into
>> PetscStackCallBLAS("LAPACKgesvd",LAPACKgrsvd_(...))
>> or should I be transposing this before passing it in?
> 
> LAPACK terminology is with respect to Fortran ordering.  There is a
> "leading dimension" parameter so that you can operate on non-contiguous
> blocks.  See KSPComputeExtremeSingularValues_GMRES for an example.
> 
>> Also for the left and right singular vector matrices that are returned by
>> this function, should I be transposing these before I interpret them as C
>> arrays?
>> 
>> I've attached my modified version of gmres.c in case this is helpful. If
>> you grep for DRL (my initials) then you'll see my changes to the code.
>> 
>> Cheers, Dave.
>> 
>> /*
>>This file implements GMRES (a Generalized Minimal Residual) method.
>>Reference:  Saad and Schultz, 1986.
>> 
>> 
>>Some comments on left vs. right preconditioning, and restarts.
>>Left and right preconditioning.
>>If right preconditioning is chosen, then the problem being solved
>>by gmres is actually
>>   My =  AB^-1 y = f
>>so the initial residual is
>>  r = f - Mx
>>Note that B^-1 y = x or y = B x, and if x is non-zero, the initial
>>residual is
>>  r = f - A x
>>The final solution is then
>>  x = B^-1 y
>> 
>>If left preconditioning is chosen, then the problem being solved is
>>   My = B^-1 A x = B^-1 f,
>>and the initial residual is
>>   r  = B^-1(f - Ax)
>> 
>>Restarts:  Restarts are basically solves with x0 not equal to zero.
>>Note that we can eliminate an extra application of B^-1 between
>>restarts as long as we don't require that the solution at the end
>>of an unsuccessful gmres iteration always be the solution x.
>> */
>> 
>> #include <../src/ksp/ksp/impls/gmres/gmresimpl.h>   /*I  "petscksp.h"  
>> I*/
>> #include  // DRL
>> #define GMRES_DELTA_DIRECTIONS 10
>> #define GMRES_DEFAULT_MAXK 30
>> static PetscErrorCode 
>> KSPGMRESUpdateHessenberg(KSP,PetscInt,PetscBool,PetscReal*);
>> static PetscErrorCode KSPGMRESBuildSoln(PetscScalar*,Vec,Vec,KSP,PetscInt);
>> 
>> PetscErrorCodeKSPSetUp_GMRES(KSP ksp)
>> {
>>  PetscInt   hh,hes,rs,cc;
>>  PetscErrorCode ierr;
>>  PetscInt   max_k,k;
>>  KSP_GMRES  *gmres = (KSP_GMRES*)ksp->data;
>> 
>>  PetscFunctionBegin;
>>  max_k = gmres->max_k;  /* restart size */
>>  hh= (max_k + 2) * (max_k + 1);
>>  hes   = (max_k + 1) * (max_k + 1);
>>  rs= (max_k + 2);
>>  cc= (max_k + 1);
>> 
>>  ierr = 
>> PetscCalloc5(hh,>hh_origin,hes,>hes_origin,rs,>rs_origin,cc,>cc_origin,cc,>ss_origin);CHKERRQ(ierr);
>>  ierr = PetscLogObjectMemory((PetscObject)ksp,(hh + hes + rs + 
>> 2*cc)*sizeof(PetscScalar));CHKERRQ(ierr);
>> 
>>  if (ksp->calc_sings) {
>>/* Allocate workspace to hold Hessenberg matrix needed by lapack */
>>ierr = PetscMalloc1((max_k + 3)*(max_k + 9),>Rsvd);CHKERRQ(ierr);
>>ierr = PetscLogObjectMemory((PetscObject)ksp,(max_k + 3)*(max_k + 
>> 9)*sizeof(PetscScalar));CHKERRQ(ierr);
>>ierr = PetscMalloc1(6*(max_k+2),>Dsvd);CHKERRQ(ierr);
>>ierr = 
>> PetscLogObjectMemory((PetscObject)ksp,6*(max_k+2)*sizeof(PetscReal));CHKERRQ(ierr);
>>  }
>> 
>>  /* Allocate array to hold pointers to user vectors.  Note that we need
>>   4 + max_k + 1 (since we need it+1 vectors, and it <= max_k) */
>>  gmres->vecs_allocated = VEC_OFFSET + 2 + max_k + gmres->nextra_vecs;
>> 
>>  ierr = PetscMalloc1(gmres->vecs_allocated,>vecs);CHKERRQ(ierr);
>>  ierr = PetscMalloc1(VEC_OFFSET+2+max_k,>user_work);CHKERRQ(ierr);
>>  ierr = PetscMalloc1(VEC_OFFSET+2+max_k,>mwork_alloc);CHKERRQ(ierr);
>>  ierr = 
>> 

Re: [petsc-users] Calling LAPACK routines from PETSc

[Jed Brown via petsc-users]

Dave Lee via petsc-users  writes:

> Hi Petsc,
>
> I'm attempting to implement a "hookstep" for the SNES trust region solver.
> Essentially what I'm trying to do is replace the solution of the least
> squares problem at the end of each GMRES solve with a modified solution
> with a norm that is constrained to be within the size of the trust region.
>
> In order to do this I need to perform an SVD on the Hessenberg matrix,
> which copying the function KSPComputeExtremeSingularValues(), I'm trying to
> do by accessing the LAPACK function dgesvd() via the PetscStackCallBLAS()
> machinery. One thing I'm confused about however is the ordering of the 2D
> arrays into and out of this function, given that that C and FORTRAN arrays
> use reverse indexing, ie: C[j+1][i+1] = F[i,j].
>
> Given that the Hessenberg matrix has k+1 rows and k columns, should I be
> still be initializing this as H[row][col] and passing this into
> PetscStackCallBLAS("LAPACKgesvd",LAPACKgrsvd_(...))
> or should I be transposing this before passing it in?

LAPACK terminology is with respect to Fortran ordering.  There is a
"leading dimension" parameter so that you can operate on non-contiguous
blocks.  See KSPComputeExtremeSingularValues_GMRES for an example.

> Also for the left and right singular vector matrices that are returned by
> this function, should I be transposing these before I interpret them as C
> arrays?
>
> I've attached my modified version of gmres.c in case this is helpful. If
> you grep for DRL (my initials) then you'll see my changes to the code.
>
> Cheers, Dave.
>
> /*
> This file implements GMRES (a Generalized Minimal Residual) method.
> Reference:  Saad and Schultz, 1986.
>
>
> Some comments on left vs. right preconditioning, and restarts.
> Left and right preconditioning.
> If right preconditioning is chosen, then the problem being solved
> by gmres is actually
>My =  AB^-1 y = f
> so the initial residual is
>   r = f - Mx
> Note that B^-1 y = x or y = B x, and if x is non-zero, the initial
> residual is
>   r = f - A x
> The final solution is then
>   x = B^-1 y
>
> If left preconditioning is chosen, then the problem being solved is
>My = B^-1 A x = B^-1 f,
> and the initial residual is
>r  = B^-1(f - Ax)
>
> Restarts:  Restarts are basically solves with x0 not equal to zero.
> Note that we can eliminate an extra application of B^-1 between
> restarts as long as we don't require that the solution at the end
> of an unsuccessful gmres iteration always be the solution x.
>  */
>
> #include <../src/ksp/ksp/impls/gmres/gmresimpl.h>   /*I  "petscksp.h"  I*/
> #include  // DRL
> #define GMRES_DELTA_DIRECTIONS 10
> #define GMRES_DEFAULT_MAXK 30
> static PetscErrorCode 
> KSPGMRESUpdateHessenberg(KSP,PetscInt,PetscBool,PetscReal*);
> static PetscErrorCode KSPGMRESBuildSoln(PetscScalar*,Vec,Vec,KSP,PetscInt);
>
> PetscErrorCodeKSPSetUp_GMRES(KSP ksp)
> {
>   PetscInt   hh,hes,rs,cc;
>   PetscErrorCode ierr;
>   PetscInt   max_k,k;
>   KSP_GMRES  *gmres = (KSP_GMRES*)ksp->data;
>
>   PetscFunctionBegin;
>   max_k = gmres->max_k;  /* restart size */
>   hh= (max_k + 2) * (max_k + 1);
>   hes   = (max_k + 1) * (max_k + 1);
>   rs= (max_k + 2);
>   cc= (max_k + 1);
>
>   ierr = 
> PetscCalloc5(hh,>hh_origin,hes,>hes_origin,rs,>rs_origin,cc,>cc_origin,cc,>ss_origin);CHKERRQ(ierr);
>   ierr = PetscLogObjectMemory((PetscObject)ksp,(hh + hes + rs + 
> 2*cc)*sizeof(PetscScalar));CHKERRQ(ierr);
>
>   if (ksp->calc_sings) {
> /* Allocate workspace to hold Hessenberg matrix needed by lapack */
> ierr = PetscMalloc1((max_k + 3)*(max_k + 9),>Rsvd);CHKERRQ(ierr);
> ierr = PetscLogObjectMemory((PetscObject)ksp,(max_k + 3)*(max_k + 
> 9)*sizeof(PetscScalar));CHKERRQ(ierr);
> ierr = PetscMalloc1(6*(max_k+2),>Dsvd);CHKERRQ(ierr);
> ierr = 
> PetscLogObjectMemory((PetscObject)ksp,6*(max_k+2)*sizeof(PetscReal));CHKERRQ(ierr);
>   }
>
>   /* Allocate array to hold pointers to user vectors.  Note that we need
>4 + max_k + 1 (since we need it+1 vectors, and it <= max_k) */
>   gmres->vecs_allocated = VEC_OFFSET + 2 + max_k + gmres->nextra_vecs;
>
>   ierr = PetscMalloc1(gmres->vecs_allocated,>vecs);CHKERRQ(ierr);
>   ierr = PetscMalloc1(VEC_OFFSET+2+max_k,>user_work);CHKERRQ(ierr);
>   ierr = PetscMalloc1(VEC_OFFSET+2+max_k,>mwork_alloc);CHKERRQ(ierr);
>   ierr = 
> PetscLogObjectMemory((PetscObject)ksp,(VEC_OFFSET+2+max_k)*(sizeof(Vec*)+sizeof(PetscInt))
>  + gmres->vecs_allocated*sizeof(Vec));CHKERRQ(ierr);
>
>   if (gmres->q_preallocate) {
> gmres->vv_allocated = VEC_OFFSET + 2 + max_k;
>
> ierr = 
> KSPCreateVecs(ksp,gmres->vv_allocated,>user_work[0],0,NULL);CHKERRQ(ierr);
> ierr = 
> PetscLogObjectParents(ksp,gmres->vv_allocated,gmres->user_work[0]);CHKERRQ(ierr);
>
> gmres->mwork_alloc[0] = gmres->vv_allocated;
> gmres->nwork_alloc= 1;
>