[petsc-users] Provide Matrix Factorization to EPS for Generalized Eigenvalue Problem

[Peetz, Darin T]

Hello,

I'm wondering how I could go about providing a matrix factorization calculated 
in Petsc to the eigenvalue routines in Slepc.  I'm trying to solve the 
eigenvalue problem for stability, where the solution to KU=F is needed to 
construct the A-matrix (K_sigma) in the eigenvalue problem.  Since the 
eigenvalue problem is generalized, it seems like the best way to solve it is to 
factorize the B-Matrix (K, same as in KU=F) with a package like MUMPS and use a 
method in Slepc such as Krylov-Schur.  Since I need to solve both KU=F and the 
eigenvalue problem, I'd like to compute the factorization of K first to solve 
KU=F, and then reuse it in the EPS routines.

I've tried using EPSGetST() and STSetKSP() to provide the KSP object that I 
used to solve KU=F, but then for some reason when I change nonzero values in K 
(but not nonzero locations) Petsc redoes the symbolic factorization when I go 
to solve KU=F again (it's part of an optimization routine, so I'm solving both 
problems, updating parts of K, and repeating).  This does provide the correct 
solution, and allows me to use the same factorization for KU=F and the 
eigenvalue problem, but the extra symbolic factorizations, while comparatively 
cheap, are unnecessary and ideally should be eliminated.  If I skip the calls 
to EPSGetST() and STSetKSP(), the symbolic factorization for the KSP object 
associated with KU=F is only performed once, as it should be.  Is there some 
option I'm overlooking, or maybe a better way to go about this?

Thanks,
Darin

Re: [petsc-users] GAMG

[Jeremy Theler]

On Mon, 2016-10-31 at 08:44 -0600, Jed Brown wrote:

> > After understanding Matt's point about the near nullspace (and reading
> > some interesting comments from Jed on scicomp stackexchange) I did built
> > my own vectors (I had to take a look at MatNullSpaceCreateRigidBody()
> > because I found out by running the code the nullspace should be an
> > orthonormal basis, it should say so in the docs).
> 
> Where?
> "vecs - the vectors that span the null space (excluding the constant 
> vector); these vectors must be orthonormal."
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatNullSpaceCreate.html

ok, I might have passed on that but I started with 

http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSetNearNullSpace.html

that says “attaches a null space to a matrix, which is often the null
space (rigid body modes) of the operator without boundary conditions
This null space will be used to provide near null space vectors to a
multigrid preconditioner built from this matrix.”

It wouldn't hurt to remind dumb users like me that “...it is often the
set of _orthonormalized_ rigid body modes...”


> And if you run in debug mode (default), as you always should until you
> are confident that your code is correct, MatNullSpaceCreate tests that
> your vectors are orthonormal.

That's how I realized I needed to normalize. Then I found
MatNullSpaceCreateRigidBody() and copied the code to orthogonalize.

Wouldn't it be better to orthonormalize inside MatSetNullSpace()? I bet
an orthogonalization from PETSc's code would beat any user-side code.

> > Now, there are some results I do not understand. I tried these six
> > combinations:
> >
> > order  near-nullspace   iterationsnorm
> > -  --   --
> > unknownexplicit 101.6e-6
> > unknownPCSetCoordinates 151.7e-7
> > unknownnone 152.4e-7
> > node   explicit fails with error -11
> > node   PCSetCoordinates fails with error -11
> > node   none 133.8e-7
> 
> Did you set a block size for the "node-based" orderings?  Are you sure
> the above is labeled correctly?  Anyway, PCSetCoordinates uses
> "node-based" ordering.  Implementation performance will generally be
> better with node-based ordering -- it has better memory streaming and
> cache behavior.

Yes. Indeed, when I save the stiffnes matrix as a binary file I get
a .info file that contains

-matload_block_size 3

The labeling is right, I re-checked. That's the funny part, I can't get
GAMG to work with PCSetCoordinates (which BTW, I think its documentation
does not address the issue of DOF ordering).

Any idea of what can be happening to me?


> The AIJ matrix format will also automatically do an "inode" optimization
> to reduce memory bandwidth and enable block smoothing (default
> configuration uses SOR smoothing).  You can use -mat_no_inode to try
> turning that off.


That option does not make any difference.

> 
> > Error -11 is 
> > PETSc's linear solver did not converge with reason
> > 'DIVERGED_PCSETUP_FAILED' (-11)
> Isn't there an actual error message?

Sorry, KSPGetConvergedReason() returns -11 and then my code prints that
error string. Find attached the output with -info.

Thanks
--
jeremy



[0] PetscInitialize(): PETSc successfully started: number of processors = 1
[0] PetscGetHostName(): Rejecting domainname, likely is NIS tom.(none)
[0] PetscInitialize(): Running on machine: tom
[0] SlepcInitialize(): SLEPc successfully started
697 3611
[0] PetscCommDuplicate(): Duplicating a communicator 2 2 max tags = 1
[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 2091 X 2091; storage space: 835506 
unneeded,74079 used
[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 81
[0] MatCheckCompressedRow(): Found the ratio (num_zerorows 0)/(num_localrows 
2091) < 0.6. Do not use CompressedRow routines.
[0] MatSeqAIJCheckInode(): Found 697 nodes of 2091. Limit used: 5. Using Inode 
routines
[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 2091 X 2091; storage space: 0 
unneeded,74079 used
[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 81
[0] MatCheckCompressedRow(): Found the ratio (num_zerorows 0)/(num_localrows 
2091) < 0.6. Do not use CompressedRow routines.
[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 2091 X 2091; storage space: 0 
unneeded,74079 used
[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 81
[0] MatCheckCompressedRow(): Found the ratio (num_zerorows 0)/(num_localrows 
2091) < 0.6. Do not use CompressedRow routines.
[0] PCSetUp(): Setting up PC for first time
[0] PCSetUp_GAMG(): level 0) N=2091, n data rows=3, n data cols=6, nnz/row 

Re: [petsc-users] GAMG

[Jed Brown]

"Kong, Fande"  writes:
> If the boundary values are not zero, no way to maintain symmetry unless we
> reduce the extra part of  the matrix. Not updating the columns is better in
> this situation.

The inhomogeneity of the boundary condition has nothing to do with operator 
symmetry.

I like this formulation for Dirichlet conditions.

https://scicomp.stackexchange.com/questions/3298/appropriate-space-for-weak-solutions-to-an-elliptical-pde-with-mixed-inhomogeneo/3300#3300


signature.asc
Description: PGP signature

Re: [petsc-users] GAMG

[Matthew Knepley]

On Mon, Oct 31, 2016 at 10:29 AM, Kong, Fande  wrote:

> On Mon, Oct 31, 2016 at 8:44 AM, Jed Brown  wrote:
>
>> Jeremy Theler  writes:
>>
>> > Hi again
>> >
>> > I have been wokring on these issues. Long story short: it is about the
>> > ordering of the unknown fields in the vector.
>> >
>> > Long story:
>> > The physics is linear elastic problem, you can see it does work with LU
>> > over a simple cube (warp the displacements to see it does represent an
>> > elastic problem, E=200e3, nu=0.3):
>> >
>> > https://caeplex.com/demo/results.php?id=5817146bdb561
>> >
>> >
>> > Say my three displacements (unknowns) are u,v,w. I can define the
>> > unknown vector as (is this called node-based ordering?)
>> >
>> > [u1 v1 w1 u2 v2 w2 ... un vn wn]^T
>> >
>> > Another option is (is this called unknown-based ordering?)
>> >
>> > [u1 u2 ... un v1 v2 ... vn w1 w2 ... wn]^T
>> >
>> >
>> > With lu/preonly the results are the same, although the stiffnes matrixes
>> > for each case are attached as PNGs. And of course, the near-nullspace
>> > vectors are different. So PCSetCoordinates() should work with one
>> > ordering and not with another one, an issue I did not take into
>> > consideration.
>> >
>> > After understanding Matt's point about the near nullspace (and reading
>> > some interesting comments from Jed on scicomp stackexchange) I did built
>> > my own vectors (I had to take a look at MatNullSpaceCreateRigidBody()
>> > because I found out by running the code the nullspace should be an
>> > orthonormal basis, it should say so in the docs).
>>
>> Where?
>>
>> "vecs   - the vectors that span the null space (excluding the constant
>> vector); these vectors must be orthonormal."
>>
>> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages
>> /Mat/MatNullSpaceCreate.html
>>
>> And if you run in debug mode (default), as you always should until you
>> are confident that your code is correct, MatNullSpaceCreate tests that
>> your vectors are orthonormal.
>>
>> > Now, there are some results I do not understand. I tried these six
>> > combinations:
>> >
>> > order  near-nullspace   iterationsnorm
>> > -  --   --
>> > unknownexplicit 101.6e-6
>> > unknownPCSetCoordinates 151.7e-7
>> > unknownnone 152.4e-7
>> > node   explicit fails with error -11
>> > node   PCSetCoordinates fails with error -11
>> > node   none 133.8e-7
>>
>> Did you set a block size for the "node-based" orderings?  Are you sure
>> the above is labeled correctly?  Anyway, PCSetCoordinates uses
>> "node-based" ordering.  Implementation performance will generally be
>> better with node-based ordering -- it has better memory streaming and
>> cache behavior.
>>
>> The AIJ matrix format will also automatically do an "inode" optimization
>> to reduce memory bandwidth and enable block smoothing (default
>> configuration uses SOR smoothing).  You can use -mat_no_inode to try
>> turning that off.
>>
>> > Error -11 is
>> > PETSc's linear solver did not converge with reason
>> > 'DIVERGED_PCSETUP_FAILED' (-11)
>>
>> Isn't there an actual error message?
>>
>> > Any explanation (for dumbs)?
>> > Another thing to take into account: I am setting the dirichlet BCs with
>> > MatZeroRows(), but I am not updating the columns to keep symmetry. Can
>> > this pose a problem for GAMG?
>>
>> Usually minor, but it is better to maintain symmetry.
>>
>
> If the boundary values are not zero, no way to maintain symmetry unless we
> reduce the extra part of  the matrix. Not updating the columns is better
> in this situation.
>

?

You just eliminate the unknowns.

   Matt


>
> Fande,
>
>
>
>


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

Re: [petsc-users] GAMG

[Kong, Fande]

On Mon, Oct 31, 2016 at 8:44 AM, Jed Brown  wrote:

> Jeremy Theler  writes:
>
> > Hi again
> >
> > I have been wokring on these issues. Long story short: it is about the
> > ordering of the unknown fields in the vector.
> >
> > Long story:
> > The physics is linear elastic problem, you can see it does work with LU
> > over a simple cube (warp the displacements to see it does represent an
> > elastic problem, E=200e3, nu=0.3):
> >
> > https://caeplex.com/demo/results.php?id=5817146bdb561
> >
> >
> > Say my three displacements (unknowns) are u,v,w. I can define the
> > unknown vector as (is this called node-based ordering?)
> >
> > [u1 v1 w1 u2 v2 w2 ... un vn wn]^T
> >
> > Another option is (is this called unknown-based ordering?)
> >
> > [u1 u2 ... un v1 v2 ... vn w1 w2 ... wn]^T
> >
> >
> > With lu/preonly the results are the same, although the stiffnes matrixes
> > for each case are attached as PNGs. And of course, the near-nullspace
> > vectors are different. So PCSetCoordinates() should work with one
> > ordering and not with another one, an issue I did not take into
> > consideration.
> >
> > After understanding Matt's point about the near nullspace (and reading
> > some interesting comments from Jed on scicomp stackexchange) I did built
> > my own vectors (I had to take a look at MatNullSpaceCreateRigidBody()
> > because I found out by running the code the nullspace should be an
> > orthonormal basis, it should say so in the docs).
>
> Where?
>
> "vecs   - the vectors that span the null space (excluding the constant
> vector); these vectors must be orthonormal."
>
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/
> MatNullSpaceCreate.html
>
> And if you run in debug mode (default), as you always should until you
> are confident that your code is correct, MatNullSpaceCreate tests that
> your vectors are orthonormal.
>
> > Now, there are some results I do not understand. I tried these six
> > combinations:
> >
> > order  near-nullspace   iterationsnorm
> > -  --   --
> > unknownexplicit 101.6e-6
> > unknownPCSetCoordinates 151.7e-7
> > unknownnone 152.4e-7
> > node   explicit fails with error -11
> > node   PCSetCoordinates fails with error -11
> > node   none 133.8e-7
>
> Did you set a block size for the "node-based" orderings?  Are you sure
> the above is labeled correctly?  Anyway, PCSetCoordinates uses
> "node-based" ordering.  Implementation performance will generally be
> better with node-based ordering -- it has better memory streaming and
> cache behavior.
>
> The AIJ matrix format will also automatically do an "inode" optimization
> to reduce memory bandwidth and enable block smoothing (default
> configuration uses SOR smoothing).  You can use -mat_no_inode to try
> turning that off.
>
> > Error -11 is
> > PETSc's linear solver did not converge with reason
> > 'DIVERGED_PCSETUP_FAILED' (-11)
>
> Isn't there an actual error message?
>
> > Any explanation (for dumbs)?
> > Another thing to take into account: I am setting the dirichlet BCs with
> > MatZeroRows(), but I am not updating the columns to keep symmetry. Can
> > this pose a problem for GAMG?
>
> Usually minor, but it is better to maintain symmetry.
>

If the boundary values are not zero, no way to maintain symmetry unless we
reduce the extra part of  the matrix. Not updating the columns is better in
this situation.

Fande,

Re: [petsc-users] --download-metis and build of stand-alone tools

[Klaij, Christiaan]

Jed,

Thanks, that line in the cmake file is exactly what I needed to
know. A petsc configure option would be nice to have, but it's
too difficult for me to do right now, I'll just hack the file
instead.

Chris


dr. ir. Christiaan Klaij  | CFD Researcher | Research & Development
MARIN | T +31 317 49 33 44 | mailto:c.kl...@marin.nl | http://www.marin.nl

MARIN news: 
http://www.marin.nl/web/News/News-items/SSSRIMARIN-seminar-November-2-Shanghai.htm


From: Jed Brown 
Sent: Monday, October 31, 2016 3:26 PM
To: Klaij, Christiaan; petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] --download-metis and build of stand-alone tools

"Klaij, Christiaan"  writes:

> Satish,
>
> I've noticed that SuperLU depends on metis and parmetis and that
> PETSc downloads the versions 5.1.0-p3 and 4.0.3-p3. These are
> different from the Karypis latest stable versions (without the
> -p3). Do I really need these -p3 versions?

They fix some portability and correctness bugs.  Those packages are
mostly unmaintained by upstream and new releases often don't fix bugs
that have reproducible test cases and patches.  So you can use the
upstream version, but it might crash due to known bugs and good luck
getting support.

> If so, after configure, compilation and installation by petsc, it
> seems that the stand-alone programs such as gpmetis are not being
> build and installed. That's a problem for me. I don't mind
> switching to the versions and config that petsc needs, but I do
> need the complete thing. Can I somehow tell petsc to also build
> the standalone tools?

PETSc only needs or wants the library.  In the pkg-metis CMakeLists.txt
file, there is a line

  #add_subdirectory("programs")

which needs to be uncommented to get the programs.  Someone could add a
conditional and plumb it into metis.py to make a PETSc configure option.
That would be a welcome contribution.

Re: [petsc-users] GAMG

[Jed Brown]

Jeremy Theler  writes:

> Hi again
>
> I have been wokring on these issues. Long story short: it is about the
> ordering of the unknown fields in the vector.
>
> Long story:
> The physics is linear elastic problem, you can see it does work with LU
> over a simple cube (warp the displacements to see it does represent an
> elastic problem, E=200e3, nu=0.3):
>
> https://caeplex.com/demo/results.php?id=5817146bdb561
>
>
> Say my three displacements (unknowns) are u,v,w. I can define the
> unknown vector as (is this called node-based ordering?)
>
> [u1 v1 w1 u2 v2 w2 ... un vn wn]^T
>
> Another option is (is this called unknown-based ordering?)
>
> [u1 u2 ... un v1 v2 ... vn w1 w2 ... wn]^T
>
>
> With lu/preonly the results are the same, although the stiffnes matrixes
> for each case are attached as PNGs. And of course, the near-nullspace
> vectors are different. So PCSetCoordinates() should work with one
> ordering and not with another one, an issue I did not take into
> consideration.
>
> After understanding Matt's point about the near nullspace (and reading
> some interesting comments from Jed on scicomp stackexchange) I did built
> my own vectors (I had to take a look at MatNullSpaceCreateRigidBody()
> because I found out by running the code the nullspace should be an
> orthonormal basis, it should say so in the docs).

Where?

"vecs   - the vectors that span the null space (excluding the constant vector); 
these vectors must be orthonormal."

https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatNullSpaceCreate.html

And if you run in debug mode (default), as you always should until you
are confident that your code is correct, MatNullSpaceCreate tests that
your vectors are orthonormal.

> Now, there are some results I do not understand. I tried these six
> combinations:
>
> order  near-nullspace   iterationsnorm
> -  --   --
> unknownexplicit 101.6e-6
> unknownPCSetCoordinates 151.7e-7
> unknownnone 152.4e-7
> node   explicit fails with error -11
> node   PCSetCoordinates fails with error -11
> node   none 133.8e-7

Did you set a block size for the "node-based" orderings?  Are you sure
the above is labeled correctly?  Anyway, PCSetCoordinates uses
"node-based" ordering.  Implementation performance will generally be
better with node-based ordering -- it has better memory streaming and
cache behavior.

The AIJ matrix format will also automatically do an "inode" optimization
to reduce memory bandwidth and enable block smoothing (default
configuration uses SOR smoothing).  You can use -mat_no_inode to try
turning that off.

> Error -11 is 
> PETSc's linear solver did not converge with reason
> 'DIVERGED_PCSETUP_FAILED' (-11)

Isn't there an actual error message?

> Any explanation (for dumbs)?
> Another thing to take into account: I am setting the dirichlet BCs with
> MatZeroRows(), but I am not updating the columns to keep symmetry. Can
> this pose a problem for GAMG?

Usually minor, but it is better to maintain symmetry.


signature.asc
Description: PGP signature

Re: [petsc-users] --download-metis and build of stand-alone tools

[Jed Brown]

"Klaij, Christiaan"  writes:

> Satish,
>
> I've noticed that SuperLU depends on metis and parmetis and that
> PETSc downloads the versions 5.1.0-p3 and 4.0.3-p3. These are
> different from the Karypis latest stable versions (without the
> -p3). Do I really need these -p3 versions?

They fix some portability and correctness bugs.  Those packages are
mostly unmaintained by upstream and new releases often don't fix bugs
that have reproducible test cases and patches.  So you can use the
upstream version, but it might crash due to known bugs and good luck
getting support.

> If so, after configure, compilation and installation by petsc, it
> seems that the stand-alone programs such as gpmetis are not being
> build and installed. That's a problem for me. I don't mind
> switching to the versions and config that petsc needs, but I do
> need the complete thing. Can I somehow tell petsc to also build
> the standalone tools?

PETSc only needs or wants the library.  In the pkg-metis CMakeLists.txt
file, there is a line

  #add_subdirectory("programs")

which needs to be uncommented to get the programs.  Someone could add a
conditional and plumb it into metis.py to make a PETSc configure option.
That would be a welcome contribution.


signature.asc
Description: PGP signature

Re: [petsc-users] --download-metis and build of stand-alone tools

[Satish Balay]

On Mon, 31 Oct 2016, Klaij, Christiaan wrote:

> Satish,
> 
> I've noticed that SuperLU depends on metis and parmetis and that
> PETSc downloads the versions 5.1.0-p3 and 4.0.3-p3. These are
> different from the Karypis latest stable versions (without the
> -p3). Do I really need these -p3 versions?

The tarballs we distribute have a bunch of patches [mostly for
protable build - but some bugfixes aswell]

https://bitbucket.org/petsc/pkg-parmetis/commits/all
https://bitbucket.org/petsc/pkg-metis/commits/all

> 
> If so, after configure, compilation and installation by petsc, it
> seems that the stand-alone programs such as gpmetis are not being
> build and installed. That's a problem for me. I don't mind
> switching to the versions and config that petsc needs, but I do
> need the complete thing. Can I somehow tell petsc to also build
> the standalone tools?

I haven't built this stuff. But presumably the process with our
tarball is similar to the one you would use with tarballs from
Karypis.

So you should be able to "cd PETSC_ARCH/externalapcakages/*metis*" and
do the build of this extra stuff as you require?

Satish


> 
> Chris
> 
> 
> dr. ir. Christiaan Klaij  | CFD Researcher | Research & Development
> MARIN | T +31 317 49 33 44 | mailto:c.kl...@marin.nl | http://www.marin.nl
> 
> MARIN news: 
> http://www.marin.nl/web/News/News-items/New-SCREENIN-JIP-open-for-participation.htm
> 
> 

[petsc-users] --download-metis and build of stand-alone tools

[Klaij, Christiaan]

Satish,

I've noticed that SuperLU depends on metis and parmetis and that
PETSc downloads the versions 5.1.0-p3 and 4.0.3-p3. These are
different from the Karypis latest stable versions (without the
-p3). Do I really need these -p3 versions?

If so, after configure, compilation and installation by petsc, it
seems that the stand-alone programs such as gpmetis are not being
build and installed. That's a problem for me. I don't mind
switching to the versions and config that petsc needs, but I do
need the complete thing. Can I somehow tell petsc to also build
the standalone tools?

Chris


dr. ir. Christiaan Klaij  | CFD Researcher | Research & Development
MARIN | T +31 317 49 33 44 | mailto:c.kl...@marin.nl | http://www.marin.nl

MARIN news: 
http://www.marin.nl/web/News/News-items/New-SCREENIN-JIP-open-for-participation.htm