Re: [petsc-users] Can we force SNES solver to do at least Newton step?

[David Knezevic via petsc-users]

>
> Hmm, it sounds like the convergence measure is bad. Maybe using a weighted
>>> norm would be better?
>>
>>
>> That's a good thought, I'd like to look into that idea too. Could you
>> please give me some guidance on how to use a weighted norm in the
>> convergence test? (Or are there any examples of doing that in the example
>> suite?)
>>
>
> I guess the idea would be the following:
>
>You have exponential terms in your equation. You get a relatively small
> residual, but unacceptable
>error. This leads me to believe that some part of the residual is
> swamping other parts you care about.
>This should be confirmed by plotting the residual for one of these
> 0-iterate cases. To prevent this
>imbalance, you would scale the residual to give more weight to the
> underrepresented parts, probably
>with an inverse of the exponential. This is just a more complicated
> form of the diagonal scaling
>commonly used to give all fields the same scale.
>

OK, I see what you mean, thanks. I'll investigate that idea. And I'll also
try SNESSetForceIteration, as suggested by Barry, since it looks like that
takes care of my initial request (to always perform at least one iteration)
without needing a custom convergence test.

Thanks,
David



On Sat, Sep 2, 2023 at 5:54 PM Matthew Knepley  wrote:
>>
>>> On Sat, Sep 2, 2023 at 5:45 PM David Knezevic <
>>> david.kneze...@akselos.com> wrote:
>>>
>>>> OK, thanks, I'll look into the custom convergence test.
>>>>
>>>> I do not understand this comment. What do you mean by "inaccurate"?
>>>>> Since we do not have the true solution, we usually say "inaccurate" for
>>>>> large residual, but you already said that the residual is small. Why would
>>>>> you want to do another iterate?
>>>>
>>>>
>>>> I agree with your comments, but the specific case I'm considering is
>>>> very numerically sensitive since it includes creep (which unfortunately has
>>>> large exponential terms in it) which is the root cause of the issues I'm
>>>> facing. Based on test cases with a known reference solution we're finding
>>>> that we get inaccurate results due to steps with "zero iterations". We can
>>>> fix this by tightening the tolerance but then we do an excessive number of
>>>> iterations in other steps. So it seems to me that ensuring that we do at
>>>> least one iteration will help here, so that's what I wanted to try.
>>>>
>>>
>>> Hmm, it sounds like the convergence measure is bad. Maybe using a
>>> weighted norm would be better?
>>>
>>>   Thanks,
>>>
>>>  Matt
>>>
>>>
>>>> Thanks again for your help.
>>>>
>>>> Best,
>>>> David
>>>>
>>>> On Sat, Sep 2, 2023 at 3:23 PM Matthew Knepley 
>>>> wrote:
>>>>
>>>>> On Sat, Sep 2, 2023 at 3:05 PM David Knezevic via petsc-users <
>>>>> petsc-users@mcs.anl.gov> wrote:
>>>>>
>>>>>> Hi all,
>>>>>>
>>>>>> I'm using the SNES solver for a plasticity model, and the issue I've
>>>>>> run into is that in some time steps the solver terminates after "NL step 
>>>>>> 0"
>>>>>> since the initial residual (based on the solution from the previous time
>>>>>> step) is below the specified tolerance.
>>>>>>
>>>>>> I gather that "NL step 0" only checks the residual and doesn't
>>>>>> actually do a Newtown update, and hence it seems that this is leading to
>>>>>> inaccurate results in some cases.
>>>>>>
>>>>>
>>>>> I do not understand this comment. What do you mean by "inaccurate"?
>>>>> Since we do not have the true solution, we usually say "inaccurate" for
>>>>> large residual, but you already said that the residual is small.
>>>>> Why would you want to do another iterate?
>>>>>
>>>>>
>>>>>> I can of course specify a smaller convergence tolerance to avoid this
>>>>>> issue, but I've found it difficult to find a smaller tolerance that works
>>>>>> well in all cases (e.g. it leads to too many iterations or
>>>>>> non-convergence). So instead what I would like to do is ensure that the
>>>>>> so

Re: [petsc-users] Can we force SNES solver to do at least Newton step?

[David Knezevic via petsc-users]

>
>
> https://petsc.org/release/manualpages/SNES/SNESSetForceIteration/#snessetforceiteration


Great, thanks!

Re: [petsc-users] Can we force SNES solver to do at least Newton step?

[David Knezevic via petsc-users]

>
> Hmm, it sounds like the convergence measure is bad. Maybe using a weighted
> norm would be better?


That's a good thought, I'd like to look into that idea too. Could you
please give me some guidance on how to use a weighted norm in the
convergence test? (Or are there any examples of doing that in the example
suite?)

Thanks,
David

On Sat, Sep 2, 2023 at 5:54 PM Matthew Knepley  wrote:

> On Sat, Sep 2, 2023 at 5:45 PM David Knezevic 
> wrote:
>
>> OK, thanks, I'll look into the custom convergence test.
>>
>> I do not understand this comment. What do you mean by "inaccurate"? Since
>>> we do not have the true solution, we usually say "inaccurate" for large
>>> residual, but you already said that the residual is small. Why would you
>>> want to do another iterate?
>>
>>
>> I agree with your comments, but the specific case I'm considering is very
>> numerically sensitive since it includes creep (which unfortunately has
>> large exponential terms in it) which is the root cause of the issues I'm
>> facing. Based on test cases with a known reference solution we're finding
>> that we get inaccurate results due to steps with "zero iterations". We can
>> fix this by tightening the tolerance but then we do an excessive number of
>> iterations in other steps. So it seems to me that ensuring that we do at
>> least one iteration will help here, so that's what I wanted to try.
>>
>
> Hmm, it sounds like the convergence measure is bad. Maybe using a weighted
> norm would be better?
>
>   Thanks,
>
>  Matt
>
>
>> Thanks again for your help.
>>
>> Best,
>> David
>>
>> On Sat, Sep 2, 2023 at 3:23 PM Matthew Knepley  wrote:
>>
>>> On Sat, Sep 2, 2023 at 3:05 PM David Knezevic via petsc-users <
>>> petsc-users@mcs.anl.gov> wrote:
>>>
>>>> Hi all,
>>>>
>>>> I'm using the SNES solver for a plasticity model, and the issue I've
>>>> run into is that in some time steps the solver terminates after "NL step 0"
>>>> since the initial residual (based on the solution from the previous time
>>>> step) is below the specified tolerance.
>>>>
>>>> I gather that "NL step 0" only checks the residual and doesn't actually
>>>> do a Newtown update, and hence it seems that this is leading to inaccurate
>>>> results in some cases.
>>>>
>>>
>>> I do not understand this comment. What do you mean by "inaccurate"?
>>> Since we do not have the true solution, we usually say "inaccurate" for
>>> large residual, but you already said that the residual is small.
>>> Why would you want to do another iterate?
>>>
>>>
>>>> I can of course specify a smaller convergence tolerance to avoid this
>>>> issue, but I've found it difficult to find a smaller tolerance that works
>>>> well in all cases (e.g. it leads to too many iterations or
>>>> non-convergence). So instead what I would like to do is ensure that the
>>>> solver does at least 1 Newton iteration instead of terminating at "NL step
>>>> 0". Is there a way to enforce this behavior, e.g. by skipping "NL step 0",
>>>> or specifying a "minimum number of iterations"? I didn't see anything like
>>>> this in the documentation, so I was wondering if there are any suggestions
>>>> on how to proceed for this.
>>>>
>>>
>>> The easiest way to do this is to write a custom convergence test that
>>> looks like this
>>>
>>> PetscErrorCode SNESConvergedDefault(SNES snes, PetscInt it, PetscReal
>>> xnorm, PetscReal snorm, PetscReal fnorm, SNESConvergedReason *reason, void
>>> *dummy)
>>> {
>>>   PetscFunctionBeginUser;
>>>   if (!it) {
>>> *reason = SNES_CONVERGED_ITERATING;
>>> PetscFunctionReturn(PETSC_SUCCESS);
>>>   }
>>>   PetscCall(SNESConvergedDefault(snes, it, xnorm, snorm, fnorm, reason,
>>> dummy));
>>>   PetscFunctionReturn(PETSC_SUCCESS);
>>> }
>>>
>>>   Thanks,
>>>
>>>  Matt
>>>
>>>
>>>> Thanks,
>>>> David
>>>>
>>>
>>>
>>> --
>>> 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/
>>> <http://www.cse.buffalo.edu/~knepley/>
>>>
>>
>
> --
> 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/
> <http://www.cse.buffalo.edu/~knepley/>
>

Re: [petsc-users] Can we force SNES solver to do at least Newton step?

[David Knezevic via petsc-users]

OK, thanks, I'll look into the custom convergence test.

I do not understand this comment. What do you mean by "inaccurate"? Since
> we do not have the true solution, we usually say "inaccurate" for large
> residual, but you already said that the residual is small. Why would you
> want to do another iterate?


I agree with your comments, but the specific case I'm considering is very
numerically sensitive since it includes creep (which unfortunately has
large exponential terms in it) which is the root cause of the issues I'm
facing. Based on test cases with a known reference solution we're finding
that we get inaccurate results due to steps with "zero iterations". We can
fix this by tightening the tolerance but then we do an excessive number of
iterations in other steps. So it seems to me that ensuring that we do at
least one iteration will help here, so that's what I wanted to try.

Thanks again for your help.

Best,
David

On Sat, Sep 2, 2023 at 3:23 PM Matthew Knepley  wrote:

> On Sat, Sep 2, 2023 at 3:05 PM David Knezevic via petsc-users <
> petsc-users@mcs.anl.gov> wrote:
>
>> Hi all,
>>
>> I'm using the SNES solver for a plasticity model, and the issue I've run
>> into is that in some time steps the solver terminates after "NL step 0"
>> since the initial residual (based on the solution from the previous time
>> step) is below the specified tolerance.
>>
>> I gather that "NL step 0" only checks the residual and doesn't actually
>> do a Newtown update, and hence it seems that this is leading to inaccurate
>> results in some cases.
>>
>
> I do not understand this comment. What do you mean by "inaccurate"? Since
> we do not have the true solution, we usually say "inaccurate" for large
> residual, but you already said that the residual is small.
> Why would you want to do another iterate?
>
>
>> I can of course specify a smaller convergence tolerance to avoid this
>> issue, but I've found it difficult to find a smaller tolerance that works
>> well in all cases (e.g. it leads to too many iterations or
>> non-convergence). So instead what I would like to do is ensure that the
>> solver does at least 1 Newton iteration instead of terminating at "NL step
>> 0". Is there a way to enforce this behavior, e.g. by skipping "NL step 0",
>> or specifying a "minimum number of iterations"? I didn't see anything like
>> this in the documentation, so I was wondering if there are any suggestions
>> on how to proceed for this.
>>
>
> The easiest way to do this is to write a custom convergence test that
> looks like this
>
> PetscErrorCode SNESConvergedDefault(SNES snes, PetscInt it, PetscReal
> xnorm, PetscReal snorm, PetscReal fnorm, SNESConvergedReason *reason, void
> *dummy)
> {
>   PetscFunctionBeginUser;
>   if (!it) {
> *reason = SNES_CONVERGED_ITERATING;
> PetscFunctionReturn(PETSC_SUCCESS);
>   }
>   PetscCall(SNESConvergedDefault(snes, it, xnorm, snorm, fnorm, reason,
> dummy));
>   PetscFunctionReturn(PETSC_SUCCESS);
> }
>
>   Thanks,
>
>  Matt
>
>
>> Thanks,
>> David
>>
>
>
> --
> 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/
> <http://www.cse.buffalo.edu/~knepley/>
>

[petsc-users] Can we force SNES solver to do at least Newton step?

[David Knezevic via petsc-users]

Hi all,

I'm using the SNES solver for a plasticity model, and the issue I've run
into is that in some time steps the solver terminates after "NL step 0"
since the initial residual (based on the solution from the previous time
step) is below the specified tolerance.

I gather that "NL step 0" only checks the residual and doesn't actually do
a Newtown update, and hence it seems that this is leading to inaccurate
results in some cases. I can of course specify a smaller convergence
tolerance to avoid this issue, but I've found it difficult to find a
smaller tolerance that works well in all cases (e.g. it leads to too many
iterations or non-convergence). So instead what I would like to do is
ensure that the solver does at least 1 Newton iteration instead of
terminating at "NL step 0". Is there a way to enforce this behavior, e.g.
by skipping "NL step 0", or specifying a "minimum number of iterations"? I
didn't see anything like this in the documentation, so I was wondering if
there are any suggestions on how to proceed for this.

Thanks,
David

Re: [petsc-users] Configure error with PETSc 3.17.4 on Ubuntu 22.04

[David Knezevic via petsc-users]

Thanks for this suggestion, using --download-mpich got it working!

Thanks,
David



On Fri, Aug 12, 2022 at 1:25 PM Satish Balay  wrote:

> We've seen issues with Ubuntu built MPI. Can you try --download-mpich and
> see if that helps?
>
> Also  you don't need --with-clanguage=cxx for complex. PETSc can be built
> with c99 compilex - i.e "--with-clanguage=c --with-scalar-type=complex"
> should work.
>
> Satish
>
> On Fri, 12 Aug 2022, David Knezevic via petsc-users wrote:
>
> > I get a configure error when building PETSc 3.17.4 on Ubuntu 22.04. The
> > configure.log is attached.
> >
> > This error seems to only occur when we use --with-clanguage=cxx. We need
> to
> > use that object when building with complex numbers.
> >
> > Guidance on how to resolve this would be most appreciated.
> >
> > Thanks,
> > David
> >
>
>

Re: [petsc-users] Configure error with PETSc 3.17.4 on Ubuntu 22.04

[David Knezevic via petsc-users]

Regarding 1.: We're using --download-mumps, so this is the mumps that was
installed during config. There was no mumps build prior to that. Does that
change anything?

Regarding 2.: We need --with-clanguage=cxx since we're also using an
external C++ library (libMesh). We've been using --with-clanguage=cxx in
PETSc with libMesh for years...



On Fri, Aug 12, 2022 at 1:24 PM Matthew Knepley  wrote:

> On Fri, Aug 12, 2022 at 1:14 PM David Knezevic via petsc-users <
> petsc-users@mcs.anl.gov> wrote:
>
>> I get a configure error when building PETSc 3.17.4 on Ubuntu 22.04. The
>> configure.log is attached.
>>
>> This error seems to only occur when we use --with-clanguage=cxx. We need
>> to use that object when building with complex numbers.
>>
>> Guidance on how to resolve this would be most appreciated.
>>
>
> 1. I think your existing MUMPS build is incompatible (but it is hard for
> us to tell). Remove it
>
>   rm -rf $PETSC_DIR/$PETSC_ARCH/externalpackages/git.mumps
>
> and reconfigure
>
> 2. You can also use C99 for complex numbers, but your application code
> might already be using C++.
>
>   Thanks,
>
>  Matt
>
>
>> Thanks,
>> David
>>
>
>
> --
> 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/
> <http://www.cse.buffalo.edu/~knepley/>
>

Re: [petsc-users] Two questions regarding SNESLinesearchPrecheck

[David Knezevic]

OK, that's clear, thank you!

David


On Tue, Mar 8, 2022 at 11:43 AM Jed Brown  wrote:

> I think SNESLineSearchApply_Basic will clarify these points. Note that the
> pre-check is applied before "taking the step", so X is x_k. You're right
> that the sign is flipped on search direction Y, as it's using -lambda below.
>
>   /* precheck */
>   ierr = SNESLineSearchPreCheck(linesearch,X,Y,_y);CHKERRQ(ierr);
>
>   /* update */
>   ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
>   [...]
>   ierr = VecCopy(W, X);CHKERRQ(ierr);
>
> David Knezevic  writes:
>
> > We're using SNESLinesearchPrecheck in order to implement some nonlinear
> > continuation methods, and we had a couple of questions:
> >
> > 1. The "search direction" that is handed to the pre-check function is
> > referred to as "y" in the documentation. We had assumed that we would
> have
> > "y = delta_x_k", where delta_x_k is as shown in the attached screenshot
> > from the PETSc manual. But after doing some testing it seems that in fact
> > we have "y = -delta_x_k", i.e. y is the NEGATIVE of the Newton step. Is
> > this correct?
> >
> > 2. In the documentation for SNESLineSearchPreCheck it says that "x" is
> the
> > "current solution". Referring again to the attached screenshot, does this
> > mean that it's x_k+1 or x_k? I assume it means x_k+1 but I wanted to
> > confirm this.
> >
> > Thanks for your help!
> >
> > Regards,
> > David
>

[petsc-users] Two questions regarding SNESLinesearchPrecheck

[David Knezevic]

We're using SNESLinesearchPrecheck in order to implement some nonlinear
continuation methods, and we had a couple of questions:

1. The "search direction" that is handed to the pre-check function is
referred to as "y" in the documentation. We had assumed that we would have
"y = delta_x_k", where delta_x_k is as shown in the attached screenshot
from the PETSc manual. But after doing some testing it seems that in fact
we have "y = -delta_x_k", i.e. y is the NEGATIVE of the Newton step. Is
this correct?

2. In the documentation for SNESLineSearchPreCheck it says that "x" is the
"current solution". Referring again to the attached screenshot, does this
mean that it's x_k+1 or x_k? I assume it means x_k+1 but I wanted to
confirm this.

Thanks for your help!

Regards,
David

Re: [petsc-users] Does MatCopy require Mats to have the same communicator?

[David Knezevic]

>
> Maybe I'm not understanding you -- matCopy does not have your
> sub-communicator so how would it create a Mat with it ...
>

MatCopy has the two Mats, so I thought I could set up the copy to be based
on the sub-communicator and MatCopy might have been set up to get the
necessary communicator from the copy Mat... that was just a guess on my
part and presumably not how it actually works...



> You probably want to use MatGetSubmatrix. This is general, but
> the output will have the same communicator, but the idle processors will be
> empty, if that is what you specify.
> Now you just need to replace the communicator with your sub communicator.
> Not sure how to do this but now you have your data in the right place at
> least.
>
> Oh, there is a method to create a Mat with a sub communicator with
> non-empty rows. Now you use this and use the communicator in the new matrix
> as your sub-communicator.
>
>
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatMPIAdjCreateNonemptySubcommMat.html
>

OK, thanks, that's helpful. Though Randall Mackie also just sent me this
link:
https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateRedundantMatrix.html

It seems like MatCreateRedundantMatrix is what I want to use, so I'll try
that out.

Best,
David





>
>
>
> On Mon, Jan 27, 2020 at 11:59 AM David Knezevic <
> david.kneze...@akselos.com> wrote:
>
>> I have a case where I'd like to copy a Mat defined on COMM_WORLD to a new
>> Mat defined on some sub-communicator. Does MatCopy support this, or would I
>> have to write a custom copy operation?
>>
>> I see here
>> <https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatConvert.html#MatConvert>
>> that MatConvert requires identical communicators, but I don't see any
>> mention of this for MatCopy, so I wanted to check.
>>
>> Thanks,
>> David
>>
>

[petsc-users] Does MatCopy require Mats to have the same communicator?

[David Knezevic]

I have a case where I'd like to copy a Mat defined on COMM_WORLD to a new
Mat defined on some sub-communicator. Does MatCopy support this, or would I
have to write a custom copy operation?

I see here

that MatConvert requires identical communicators, but I don't see any
mention of this for MatCopy, so I wanted to check.

Thanks,
David

Re: [petsc-users] Is it possible to update SNES tolerances?

[David Knezevic via petsc-users]

On Mon, Jun 10, 2019 at 2:25 PM Jed Brown  wrote:

> David Knezevic via petsc-users  writes:
>
> > I'm doing load stepping with SNES, where I do a SNES solve for each load
> > step. Ideally the convergence tolerances for each load step would be set
> > based on the norm of the load in the current step. As a result I would
> like
> > to be able to update the absolute and relative tolerances on a SNES.
> >
> > I tried to call SNESSetTolerances at the beginning of each load step, but
> > it seems (based on the -snes_view output) that the tolerances do not get
> > updated. Is there a way to update the tolerances?
>
> Did you also have SNESSetFromOptions and specify tolerances on the
> command line?  SNESSetTolerances should be used unless it gets
> overwritten later by a process such as SNESSetFromOptions.
>

Yes, you're right! Oops, I was overwriting the value... all good, sorry for
the spam.

Best regards,
David

[petsc-users] Is it possible to update SNES tolerances?

[David Knezevic via petsc-users]

I'm doing load stepping with SNES, where I do a SNES solve for each load
step. Ideally the convergence tolerances for each load step would be set
based on the norm of the load in the current step. As a result I would like
to be able to update the absolute and relative tolerances on a SNES.

I tried to call SNESSetTolerances at the beginning of each load step, but
it seems (based on the -snes_view output) that the tolerances do not get
updated. Is there a way to update the tolerances?

Thank you,
David

Re: [petsc-users] Failure of MUMPS

[David Knezevic]

On Fri, Oct 5, 2018 at 9:16 PM Mike Wick 
wrote:

> Hi Hong:
>
> There is no explicit error message actually. The solver converge in 0
> iteration and returns an nan number. Looks like a zero pivoting problem.
>
> Mike
>

I use MatMumpsGetInfo

to
get error info returned by MUMPS, e.g.

PetscInt info_1;
MatMumpsGetInfo(pc_mat, 1, _1);

Then you can check the value of info_1 to get error diagnostics. One error
that I run into sometimes is when info_1 is -9, in which case I increase
icntl_14 and try again.

Best,
David



>
> On Fri, Oct 5, 2018 at 6:12 PM Zhang, Hong  wrote:
>
>> Mike:
>>
>>> Hello PETSc team:
>>>
>>> I am trying to solve a PDE problem with high-order finite elements. The
>>> matrix is getting denser and my experience is that MUMPS just outperforms
>>> iterative solvers.
>>>
>>> For certain problems, MUMPS just fail in the middle for no clear reason.
>>> I just wander if there is any suggestion to improve the robustness of
>>> MUMPS? Or in general, any suggestion for interative solver with very
>>> high-order finite elements?
>>>
>>
>> What error message do you get when MUMPS fails? Out of memory, zero
>> pivoting, or something?
>>  Hong
>>
>

Re: [petsc-users] Changing the jacobian matrix during SNES iterations

[David Knezevic]

On Wed, Jul 18, 2018 at 3:34 PM, Matthew Knepley  wrote:

> On Wed, Jul 18, 2018 at 3:25 PM David Knezevic 
> wrote:
>
>> On Wed, Jul 18, 2018 at 1:59 PM, Matthew Knepley 
>> wrote:
>>
>>> On Wed, Jul 18, 2018 at 1:31 PM David Knezevic <
>>> david.kneze...@akselos.com> wrote:
>>>
>>>> I'm using SNES for a finite element contact solve, in which the
>>>> sparsity pattern of the jacobian can change from one Newton iteration to
>>>> the next (since the nodes on the contact surface move).
>>>>
>>>> In order to handle this I figured the best way would be to destroy the
>>>> jacobian matrix and re-allocate it with a new sparsity pattern inside each
>>>> call to FormJacobian, does that seem like a reasonable approach in this
>>>> context?
>>>>
>>>
>>> Yes.
>>>
>>>
>>>> Also, I recall from an earlier discussion that this matrix
>>>> re-allocation inside FormJacobian is supported by SNES, but I just wanted
>>>> to confirm that?
>>>>
>>>
>>> Yes.
>>>
>>>
>>>> Also, I was wondering if there is any example where the matrix is
>>>> re-allocated inside SNES iterations so that I can make sure that I do it
>>>> correctly?
>>>>
>>>
>>> No, unfortunately. Contributions always welcome :)
>>>
>>
>>
>> OK, as a test case I'd like to modify snes/tutorials/ex1.c to destroy and
>> reallocate the jacobian matrix with the same sparsity pattern
>> inside FormJacobian1 (once I can do that with the same sparsity pattern,
>> then it should be straightforward to do the same thing with a modified
>> sparsity pattern). To do that, I tried adding the following code inside
>> FormJacobian1:
>>
>>   ierr = MatDestroy();
>>
>
> You do not want to destroy the matrix. There would be no way to get
> another Mat back out of the function.
>
>
>>   ierr = MatCreate(PETSC_COMM_WORLD,);CHKERRQ(ierr);
>>
>
> And do not recreate it.
>
>
>>   ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr);
>>
>
> Nor reset the sizes. However, you do want to call
>
>   MatSetPreallocationXAIJ();
>
> which should work
>
>
>>   ierr = MatSetFromOptions(B);CHKERRQ(ierr);
>>   ierr = MatSetUp(B);CHKERRQ(ierr);
>>
>
>   THanks,
>
> Matt
>

Thanks! Calling MatSetPreallocationXAIJ inside FormJacobian works for me.

So I guess this means that we create a new sparsity pattern "from scratch" each
time we call MatSetPreallocationXAIJ and set the values in the matrix?

Thanks,
David

Re: [petsc-users] Changing the jacobian matrix during SNES iterations

[David Knezevic]

On Wed, Jul 18, 2018 at 1:59 PM, Matthew Knepley  wrote:

> On Wed, Jul 18, 2018 at 1:31 PM David Knezevic 
> wrote:
>
>> I'm using SNES for a finite element contact solve, in which the sparsity
>> pattern of the jacobian can change from one Newton iteration to the next
>> (since the nodes on the contact surface move).
>>
>> In order to handle this I figured the best way would be to destroy the
>> jacobian matrix and re-allocate it with a new sparsity pattern inside each
>> call to FormJacobian, does that seem like a reasonable approach in this
>> context?
>>
>
> Yes.
>
>
>> Also, I recall from an earlier discussion that this matrix re-allocation
>> inside FormJacobian is supported by SNES, but I just wanted to confirm that?
>>
>
> Yes.
>
>
>> Also, I was wondering if there is any example where the matrix is
>> re-allocated inside SNES iterations so that I can make sure that I do it
>> correctly?
>>
>
> No, unfortunately. Contributions always welcome :)
>


OK, as a test case I'd like to modify snes/tutorials/ex1.c to destroy and
reallocate the jacobian matrix with the same sparsity pattern
inside FormJacobian1 (once I can do that with the same sparsity pattern,
then it should be straightforward to do the same thing with a modified
sparsity pattern). To do that, I tried adding the following code inside
FormJacobian1:

  ierr = MatDestroy();
  ierr = MatCreate(PETSC_COMM_WORLD,);CHKERRQ(ierr);
  ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr);
  ierr = MatSetFromOptions(B);CHKERRQ(ierr);
  ierr = MatSetUp(B);CHKERRQ(ierr);

That gives a segfault, so I gather something else is needed here, any
suggestions on what I need to do? I guess I need to do something else in
order to make sure that the SNES uses the newly created matrix (e.g. update
the original Mat J)?

Thanks,
David

[petsc-users] Changing the jacobian matrix during SNES iterations

[David Knezevic]

I'm using SNES for a finite element contact solve, in which the sparsity
pattern of the jacobian can change from one Newton iteration to the next
(since the nodes on the contact surface move).

In order to handle this I figured the best way would be to destroy the
jacobian matrix and re-allocate it with a new sparsity pattern inside each
call to FormJacobian, does that seem like a reasonable approach in this
context? Also, I recall from an earlier discussion that this matrix
re-allocation inside FormJacobian is supported by SNES, but I just wanted
to confirm that?

Also, I was wondering if there is any example where the matrix is
re-allocated inside SNES iterations so that I can make sure that I do it
correctly?

Best regards,
David

Re: [petsc-users] Replicating a hang with MUMPS

[David Knezevic]

On Wed, Jun 27, 2018 at 3:12 PM, Smith, Barry F.  wrote:

>
>   David,
>
>  This is ugly but should work. BEFORE reading in the matrix and right
> hand side set the LOCAL sizes for the matrix and vector. This way you can
> control exactly which rows go on which process. Note you will have to have
> your own mechanism to know what the local sizes should be (for example have
> the original program print out the sizes and just cut and paste them into
> your copy of ex10.c) PETSc doesn't provide an automatic way to do this (nor
> should it).
>
> Barry
>


Thanks Barry and Stefano, the approach you both suggested (calling
MatSetSizes and VecSetSizes before reading in) was what I was looking for.

Best,
David




> On Jun 27, 2018, at 1:36 PM, David Knezevic 
> wrote:
> >
> > I ran into a case where using MUMPS (called via "-ksp_type preonly
> -pc_type lu -pc_factor_mat_solver_package mumps") for a particular solve
> hangs indefinitely with 24 MPI processes (but it works fine with other
> numbers of processes). The stack trace when killing the job is below, in
> case that gives any clue as to what is wrong.
> >
> > I'm trying to replicate this with a simple test case. I wrote out the
> matrix and right-hand side to disk using MatView and VecView, and then I
> modified ksp ex10 to read in these files and solve with 24 cores. However,
> that did not replicate the error, so I think I also need to make sure that
> I use the same number of rows per process in the test case as in the case
> that hung. As a result I'm wondering if there is a way to modify the
> parallel layout of the matrix and vector after I read them in?
> >
> > Also, if there are any other suggestions about reproducing or debugging
> this issue, please let me know!
> >
> > Best,
> > David
> >
> > 
> >
> > #0  0x7fb12bf0e74d in poll () at ../sysdeps/unix/syscall-templa
> te.S:84
> > #1  0x7fb126262e58 in ?? () from /usr/lib/libopen-pal.so.13
> > #2  0x7fb1262596fb in opal_libevent2021_event_base_loop () from
> /usr/lib/libopen-pal.so.13
> > #3  0x7fb126223238 in opal_progress () from
> /usr/lib/libopen-pal.so.13
> > #4  0x7fb12cef53db in ompi_request_default_test () from
> /usr/lib/libmpi.so.12
> > #5  0x7fb12cf21d61 in PMPI_Test () from /usr/lib/libmpi.so.12
> > #6  0x7fb127a5b939 in pmpi_test__ () from /usr/lib/libmpi_mpifh.so.12
> > #7  0x7fb132888d87 in dmumps_try_recvtreat (comm_load=8,
> ass_irecv=40, blocking=.FALSE., set_irecv=.TRUE., message_received=.FALSE.,
> msgsou=-1, msgtag=-1, status=..., bufr=..., lbufr=401408,
> lbufr_bytes=1605629, procnode_steps=..., posfac=410095, iwpos=3151,
> iwposcb=30557,
> > iptrlu=1536548, lrlu=1126454, lrlus=2864100, n=30675, iw=...,
> liw=39935, a=..., la=3367108, ptrist=..., ptlust=..., ptrfac=...,
> ptrast=..., step=..., pimaster=..., pamaster=..., nstk_s=..., comp=0,
> iflag=0, ierror=0, comm=7, nbprocfils=..., ipool=..., lpool=48, leaf=2,
> > nbfin=90, myid=33, slavef=90, root=..., opassw=353031,
> opeliw=700399235, itloc=..., rhs_mumps=..., fils=..., ptrarw=...,
> ptraiw=..., intarr=..., dblarr=..., icntl=..., keep=..., keep8=...,
> dkeep=..., nd=..., frere=..., lptrar=30675, nelt=1, frtptr=..., frtelt=...,
> > istep_to_iniv2=..., tab_pos_in_pere=...,
> stack_right_authorized=.TRUE., lrgroups=...) at dfac_process_message.F:646
> > #8  0x7fb1328cfcd1 in dmumps_fac_par_m::dmumps_fac_par (n=30675,
> iw=..., liw=39935, a=..., la=3367108, nstk_steps=..., nbprocfils=...,
> nd=..., fils=..., step=..., frere=..., dad=..., cand=...,
> istep_to_iniv2=..., tab_pos_in_pere=..., maxfrt=0, ntotpv=0, nmaxnpiv=150,
> > ptrist=..., ptrast=..., pimaster=..., pamaster=..., ptrarw=...,
> ptraiw=..., itloc=..., rhs_mumps=..., ipool=..., lpool=48, rinfo=...,
> posfac=410095, iwpos=3151, lrlu=1126454, iptrlu=1536548, lrlus=2864100,
> leaf=2, nbroot=1, nbrtot=90, uu=0.01, icntl=..., ptlust=..., ptrfac=...,
> > nsteps=1, info=..., keep=..., keep8=..., procnode_steps=...,
> slavef=90, myid=33, comm_nodes=7, myid_nodes=33, bufr=..., lbufr=401408,
> lbufr_bytes=1605629, intarr=..., dblarr=..., root=..., perm=..., nelt=1,
> frtptr=..., frtelt=..., lptrar=30675, comm_load=8, ass_irecv=40,
> > seuil=0, seuil_ldlt_niv2=0, mem_distrib=..., ne=..., dkeep=...,
> pivnul_list=..., lpn_list=1, lrgroups=...) at dfac_par_m.F:207
> > #9  0x7fb13287f875 in dmumps_fac_b (n=30675, nsteps=1, a=...,
> la=3367108, iw=..., liw=39935, sym_perm=..., na=..., lna=47, ne_steps=...,
> nfsiz=..., fils=..., step=..., frere=..., dad=..., cand=...,
> istep_to_iniv2=..., tab_pos_in_pere=..., ptrar=..., ldptrar=30675,
> 

[petsc-users] Replicating a hang with MUMPS

[David Knezevic]

I ran into a case where using MUMPS (called via "-ksp_type preonly -pc_type
lu -pc_factor_mat_solver_package mumps") for a particular solve hangs
indefinitely with 24 MPI processes (but it works fine with other numbers of
processes). The stack trace when killing the job is below, in case that
gives any clue as to what is wrong.

I'm trying to replicate this with a simple test case. I wrote out the
matrix and right-hand side to disk using MatView and VecView, and then I
modified ksp ex10 to read in these files and solve with 24 cores. However,
that did not replicate the error, so I think I also need to make sure that
I use the same number of rows per process in the test case as in the case
that hung. As a result I'm wondering if there is a way to modify the
parallel layout of the matrix and vector after I read them in?

Also, if there are any other suggestions about reproducing or debugging
this issue, please let me know!

Best,
David



#0  0x7fb12bf0e74d in poll () at ../sysdeps/unix/syscall-template.S:84
#1  0x7fb126262e58 in ?? () from /usr/lib/libopen-pal.so.13
#2  0x7fb1262596fb in opal_libevent2021_event_base_loop () from
/usr/lib/libopen-pal.so.13
#3  0x7fb126223238 in opal_progress () from /usr/lib/libopen-pal.so.13
#4  0x7fb12cef53db in ompi_request_default_test () from
/usr/lib/libmpi.so.12
#5  0x7fb12cf21d61 in PMPI_Test () from /usr/lib/libmpi.so.12
#6  0x7fb127a5b939 in pmpi_test__ () from /usr/lib/libmpi_mpifh.so.12
#7  0x7fb132888d87 in dmumps_try_recvtreat (comm_load=8, ass_irecv=40,
blocking=.FALSE., set_irecv=.TRUE., message_received=.FALSE., msgsou=-1,
msgtag=-1, status=..., bufr=..., lbufr=401408, lbufr_bytes=1605629,
procnode_steps=..., posfac=410095, iwpos=3151, iwposcb=30557,
iptrlu=1536548, lrlu=1126454, lrlus=2864100, n=30675, iw=...,
liw=39935, a=..., la=3367108, ptrist=..., ptlust=..., ptrfac=...,
ptrast=..., step=..., pimaster=..., pamaster=..., nstk_s=..., comp=0,
iflag=0, ierror=0, comm=7, nbprocfils=..., ipool=..., lpool=48, leaf=2,
nbfin=90, myid=33, slavef=90, root=..., opassw=353031,
opeliw=700399235, itloc=..., rhs_mumps=..., fils=..., ptrarw=...,
ptraiw=..., intarr=..., dblarr=..., icntl=..., keep=..., keep8=...,
dkeep=..., nd=..., frere=..., lptrar=30675, nelt=1, frtptr=..., frtelt=...,
istep_to_iniv2=..., tab_pos_in_pere=..., stack_right_authorized=.TRUE.,
lrgroups=...) at dfac_process_message.F:646
#8  0x7fb1328cfcd1 in dmumps_fac_par_m::dmumps_fac_par (n=30675,
iw=..., liw=39935, a=..., la=3367108, nstk_steps=..., nbprocfils=...,
nd=..., fils=..., step=..., frere=..., dad=..., cand=...,
istep_to_iniv2=..., tab_pos_in_pere=..., maxfrt=0, ntotpv=0, nmaxnpiv=150,
ptrist=..., ptrast=..., pimaster=..., pamaster=..., ptrarw=...,
ptraiw=..., itloc=..., rhs_mumps=..., ipool=..., lpool=48, rinfo=...,
posfac=410095, iwpos=3151, lrlu=1126454, iptrlu=1536548, lrlus=2864100,
leaf=2, nbroot=1, nbrtot=90, uu=0.01, icntl=..., ptlust=..., ptrfac=...,
nsteps=1, info=..., keep=..., keep8=..., procnode_steps=..., slavef=90,
myid=33, comm_nodes=7, myid_nodes=33, bufr=..., lbufr=401408,
lbufr_bytes=1605629, intarr=..., dblarr=..., root=..., perm=..., nelt=1,
frtptr=..., frtelt=..., lptrar=30675, comm_load=8, ass_irecv=40,
seuil=0, seuil_ldlt_niv2=0, mem_distrib=..., ne=..., dkeep=...,
pivnul_list=..., lpn_list=1, lrgroups=...) at dfac_par_m.F:207
#9  0x7fb13287f875 in dmumps_fac_b (n=30675, nsteps=1, a=...,
la=3367108, iw=..., liw=39935, sym_perm=..., na=..., lna=47, ne_steps=...,
nfsiz=..., fils=..., step=..., frere=..., dad=..., cand=...,
istep_to_iniv2=..., tab_pos_in_pere=..., ptrar=..., ldptrar=30675,
ptrist=...,
ptlust_s=..., ptrfac=..., iw1=..., iw2=..., itloc=..., rhs_mumps=...,
pool=..., lpool=48, cntl1=0.01, icntl=..., info=..., rinfo=..., keep=...,
keep8=..., procnode_steps=..., slavef=90, comm_nodes=7, myid=33,
myid_nodes=33, bufr=..., lbufr=401408, lbufr_bytes=1605629,
intarr=..., dblarr=..., root=..., nelt=1, frtptr=..., frtelt=...,
comm_load=8, ass_irecv=40, seuil=0, seuil_ldlt_niv2=0, mem_distrib=...,
dkeep=..., pivnul_list=..., lpn_list=1, lrgroups=...) at dfac_b.F:167
#10 0x7fb1328419ed in dmumps_fac_driver (id=) at
dfac_driver.F:2291
#11 0x7fb1327ff6dc in dmumps (id=) at
dmumps_driver.F:1686
#12 0x7fb1327faf0a in dmumps_f77 (job=2, sym=0, par=1, comm_f77=5,
n=30675, icntl=..., cntl=..., keep=..., dkeep=..., keep8=..., nz=0, nnz=0,
irn=..., irnhere=0, jcn=..., jcnhere=0, a=..., ahere=0, nz_loc=622296,
nnz_loc=0, irn_loc=..., irn_lochere=1, jcn_loc=...,
jcn_lochere=1, a_loc=..., a_lochere=1, nelt=0, eltptr=...,
eltptrhere=0, eltvar=..., eltvarhere=0, a_elt=..., a_elthere=0,
perm_in=..., perm_inhere=0, rhs=..., rhshere=0, redrhs=..., redrhshere=0,
info=..., rinfo=..., infog=..., rinfog=..., deficiency=0, lwk_user=0,
size_schur=0, listvar_schur=..., listvar_schurhere=0, schur=...,
schurhere=0, wk_user=..., wk_userhere=0, colsca=..., colscahere=0,

Re: [petsc-users] DIVERGED_DTOL with SNES after switching to 3.8

[David Knezevic]

That works, thanks!

David

On Sat, Oct 28, 2017 at 10:44 AM, Smith, Barry F. <bsm...@mcs.anl.gov>
wrote:

>
>   Run with something like   -snes_divergence_tolerance 1e20
>
>   Barry
>
>
>
>
>
>
>
> > On Oct 28, 2017, at 9:31 AM, David Knezevic <david.kneze...@akselos.com>
> wrote:
> >
> > Hi all,
> >
> > I updated to the head of the maint branch (b5b0337). Previously (with
> version 3.7), I had a SNES test case that converged like this:
> >
> >   NL step  0, |residual|_2 = 9.074962e-03
> >   NL step  1, |residual|_2 = 5.516783e+03
> >   NL step  2, |residual|_2 = 1.370884e-12
> > Number of nonlinear iterations: 2
> > Nonlinear solver convergence/divergence reason: CONVERGED_FNORM_RELATIVE
> >
> > and now I get the following behavior:
> >
> >   NL step  0, |residual|_2 = 9.074962e-03
> >   NL step  1, |residual|_2 = 5.516783e+03
> > Number of nonlinear iterations: 1
> > Nonlinear solver convergence/divergence reason: DIVERGED_DTOL
> >
> > I didn't change anything in my code, so I guess the default DTOL options
> in SNES have changed? If someone could suggest how I can recover the old
> behavior, that'd be great.
> >
> > (Note that this is a contact problem using an augmented Lagrangian
> penalty method, so the convergence is pretty nasty and it's hard to avoid
> these large jumps in the residual, so I don't want to get the DIVERGED_DTOL
> error in this case.)
> >
> > Thanks,
> > David
>
>

[petsc-users] DIVERGED_DTOL with SNES after switching to 3.8

[David Knezevic]

Hi all,

I updated to the head of the maint branch (b5b0337). Previously (with
version 3.7), I had a SNES test case that converged like this:

  NL step  0, |residual|_2 = 9.074962e-03
  NL step  1, |residual|_2 = 5.516783e+03
  NL step  2, |residual|_2 = 1.370884e-12
Number of nonlinear iterations: 2
Nonlinear solver convergence/divergence reason: CONVERGED_FNORM_RELATIVE

and now I get the following behavior:

  NL step  0, |residual|_2 = 9.074962e-03
  NL step  1, |residual|_2 = 5.516783e+03
Number of nonlinear iterations: 1
Nonlinear solver convergence/divergence reason: DIVERGED_DTOL

I didn't change anything in my code, so I guess the default DTOL options in
SNES have changed? If someone could suggest how I can recover the old
behavior, that'd be great.

(Note that this is a contact problem using an augmented Lagrangian penalty
method, so the convergence is pretty nasty and it's hard to avoid these
large jumps in the residual, so I don't want to get the DIVERGED_DTOL error
in this case.)

Thanks,
David

Re: [petsc-users] Understanding inner vs outer fieldsplit convergence

[David Knezevic]

 type: seqsbaij
rows=30, cols=30
package used to perform factorization: petsc
total: nonzeros=285, allocated nonzeros=285
total number of mallocs used during MatSetValues
calls =0
block size is 1
linear system matrix = precond matrix:
Mat Object:(fieldsplit_block_1_)
  1 MPI processes
  type: seqaij
  rows=30, cols=30
  total: nonzeros=468, allocated nonzeros=468
  total number of mallocs used during MatSetValues calls =0
using I-node routines: found 10 nodes, limit used is 5
A01
  Mat Object:   1 MPI processes
type: seqaij
rows=30, cols=10512
total: nonzeros=1872, allocated nonzeros=1872
total number of mallocs used during MatSetValues calls =0
  using I-node routines: found 9 nodes, limit used is 5
Mat Object:(fieldsplit_block_2_) 1 MPI processes
  type: seqaij
  rows=10512, cols=10512
  total: nonzeros=372636, allocated nonzeros=372636
  total number of mallocs used during MatSetValues calls =0
using I-node routines: found 3504 nodes, limit used is 5
  linear system matrix = precond matrix:
  Mat Object:  ()   1 MPI processes
type: seqaij
rows=10542, cols=10542
total: nonzeros=376848, allocated nonzeros=377064
total number of mallocs used during MatSetValues calls =0
  using I-node routines: found 3514 nodes, limit used is 5





>
>
>
> > On Jan 26, 2017, at 1:26 PM, David Knezevic <david.kneze...@akselos.com>
> wrote:
> >
> > I'm exploring fieldsplit with Schur (this continues from some emails I
> sent a few weeks ago about this topic), and I had a quick question about
> the inner vs outer convergence.
> >
> > I've pasted the output below from "-ksp_monitor 
> > -fieldsplit_FE_split_ksp_monitor",
> and I'm just wondering about why the second outer iteration has two inner
> iteration loops, whereas all the other outer iterations have one inner
> iteration loop? I assume it is something to do with a convergence
> tolerance, but it's not clear to me which tolerance would control that.
> >
> > Thanks,
> > David
> >
> > --
> >
> > Residual norms for fieldsplit_FE_split_ solve.
> > 0 KSP Residual norm 4.742303891408e+01
> > 1 KSP Residual norm 2.909253505630e-01
> > 2 KSP Residual norm 9.891933795059e-02
> > 3 KSP Residual norm 7.147789520745e-02
> > 4 KSP Residual norm 1.668752967907e-02
> > 5 KSP Residual norm 5.019869896662e-03
> > 6 KSP Residual norm 2.848579237244e-03
> > 7 KSP Residual norm 2.847897269641e-03
> > 8 KSP Residual norm 2.840502392022e-03
> > 9 KSP Residual norm 2.831875522381e-03
> >10 KSP Residual norm 2.688309287993e-03
> >11 KSP Residual norm 1.351494303229e-03
> >12 KSP Residual norm 1.350874246297e-03
> >13 KSP Residual norm 9.154691604943e-06
> >   0 KSP Residual norm 2.254632353893e+02
> > Residual norms for fieldsplit_FE_split_ solve.
> > 0 KSP Residual norm 4.742303891408e+01
> > 1 KSP Residual norm 2.909253505630e-01
> > 2 KSP Residual norm 9.891933795059e-02
> > 3 KSP Residual norm 7.147789520745e-02
> > 4 KSP Residual norm 1.668752967907e-02
> > 5 KSP Residual norm 5.019869896662e-03
> > 6 KSP Residual norm 2.848579237244e-03
> > 7 KSP Residual norm 2.847897269641e-03
> > 8 KSP Residual norm 2.840502392022e-03
> > 9 KSP Residual norm 2.831875522381e-03
> >10 KSP Residual norm 2.688309287993e-03
> >11 KSP Residual norm 1.351494303229e-03
> >12 KSP Residual norm 1.350874246297e-03
> >13 KSP Residual norm 9.154691604943e-06
> > Residual norms for fieldsplit_FE_split_ solve.
> > 0 KSP Residual norm 1.554697370480e-05
> > 1 KSP Residual norm 1.554471967929e-05
> > 2 KSP Residual norm 1.551293889691e-05
> > 3 KSP Residual norm 8.031337431574e-06
> > 4 KSP Residual norm 4.137185786243e-06
> > 5 KSP Residual norm 4.066606123330e-06
> > 6 KSP Residual norm 4.051107282928e-06
> > 7 KSP Residual norm 4.047442850256e-06
> > 8 KSP Residual norm 4.047129984657e-06
> > 9 KSP Residual norm 4.030697964677e-06
> >10 KSP Residual norm 2.882383190940e-06
> >11 KSP Residual norm 3.325005138484e-07
> >12 KSP Residual norm 

[petsc-users] Understanding inner vs outer fieldsplit convergence

[David Knezevic]

I'm exploring fieldsplit with Schur (this continues from some emails I sent
a few weeks ago about this topic), and I had a quick question about the
inner vs outer convergence.

I've pasted the output below from "-ksp_monitor
-fieldsplit_FE_split_ksp_monitor", and I'm just wondering about why the
second outer iteration has two inner iteration loops, whereas all the other
outer iterations have one inner iteration loop? I assume it is something to
do with a convergence tolerance, but it's not clear to me which tolerance
would control that.

Thanks,
David

--

Residual norms for fieldsplit_FE_split_ solve.
0 KSP Residual norm 4.742303891408e+01
1 KSP Residual norm 2.909253505630e-01
2 KSP Residual norm 9.891933795059e-02
3 KSP Residual norm 7.147789520745e-02
4 KSP Residual norm 1.668752967907e-02
5 KSP Residual norm 5.019869896662e-03
6 KSP Residual norm 2.848579237244e-03
7 KSP Residual norm 2.847897269641e-03
8 KSP Residual norm 2.840502392022e-03
9 KSP Residual norm 2.831875522381e-03
   10 KSP Residual norm 2.688309287993e-03
   11 KSP Residual norm 1.351494303229e-03
   12 KSP Residual norm 1.350874246297e-03
   13 KSP Residual norm 9.154691604943e-06
  0 KSP Residual norm 2.254632353893e+02
Residual norms for fieldsplit_FE_split_ solve.
0 KSP Residual norm 4.742303891408e+01
1 KSP Residual norm 2.909253505630e-01
2 KSP Residual norm 9.891933795059e-02
3 KSP Residual norm 7.147789520745e-02
4 KSP Residual norm 1.668752967907e-02
5 KSP Residual norm 5.019869896662e-03
6 KSP Residual norm 2.848579237244e-03
7 KSP Residual norm 2.847897269641e-03
8 KSP Residual norm 2.840502392022e-03
9 KSP Residual norm 2.831875522381e-03
   10 KSP Residual norm 2.688309287993e-03
   11 KSP Residual norm 1.351494303229e-03
   12 KSP Residual norm 1.350874246297e-03
   13 KSP Residual norm 9.154691604943e-06
Residual norms for fieldsplit_FE_split_ solve.
0 KSP Residual norm 1.554697370480e-05
1 KSP Residual norm 1.554471967929e-05
2 KSP Residual norm 1.551293889691e-05
3 KSP Residual norm 8.031337431574e-06
4 KSP Residual norm 4.137185786243e-06
5 KSP Residual norm 4.066606123330e-06
6 KSP Residual norm 4.051107282928e-06
7 KSP Residual norm 4.047442850256e-06
8 KSP Residual norm 4.047129984657e-06
9 KSP Residual norm 4.030697964677e-06
   10 KSP Residual norm 2.882383190940e-06
   11 KSP Residual norm 3.325005138484e-07
   12 KSP Residual norm 2.107354774516e-07
   13 KSP Residual norm 2.107005548204e-07
   14 KSP Residual norm 4.399320792736e-08
   15 KSP Residual norm 4.236902403786e-08
   16 KSP Residual norm 2.932877082709e-08
   17 KSP Residual norm 3.881909203171e-09
   18 KSP Residual norm 1.107791399514e-09
   19 KSP Residual norm 2.645048006100e-11
  1 KSP Residual norm 8.266776463696e-01
Residual norms for fieldsplit_FE_split_ solve.
0 KSP Residual norm 9.262528453386e-08
1 KSP Residual norm 5.683232925010e-10
2 KSP Residual norm 1.915223168286e-10
3 KSP Residual norm 1.397893184942e-10
4 KSP Residual norm 1.691441435404e-11
5 KSP Residual norm 6.138315243419e-12
6 KSP Residual norm 5.576043830003e-12
7 KSP Residual norm 5.574440028225e-12
8 KSP Residual norm 5.559544964428e-12
9 KSP Residual norm 5.539862581746e-12
   10 KSP Residual norm 5.258329460152e-12
   11 KSP Residual norm 2.643581511791e-12
   12 KSP Residual norm 2.641293392449e-12
   13 KSP Residual norm 2.354608977643e-14
  2 KSP Residual norm 4.450925351013e-07
Residual norms for fieldsplit_FE_split_ solve.
0 KSP Residual norm 6.653681330477e-14
1 KSP Residual norm 6.650750698147e-14
2 KSP Residual norm 6.23464526e-14
3 KSP Residual norm 2.026817941567e-14
4 KSP Residual norm 9.604999144183e-15
5 KSP Residual norm 9.208296307424e-15
6 KSP Residual norm 9.196769686859e-15
7 KSP Residual norm 9.185058975459e-15
8 KSP Residual norm 9.180207477303e-15
9 KSP Residual norm 8.991574890909e-15
   10 KSP Residual norm 8.032736869820e-15
   11 KSP Residual norm 1.536409278928e-15
   12 KSP Residual norm 1.177374264280e-15
   13 KSP Residual norm 1.175712092044e-15
   14 KSP Residual norm 2.572275406087e-16
   15 KSP Residual norm 2.548423809711e-16
   16 KSP Residual norm 8.616505207588e-17
   17 KSP Residual norm 7.563053994201e-18
   18 KSP Residual norm 6.807636198601e-18
   19 KSP Residual norm 9.747028518744e-19
   20 KSP Residual norm 2.419807103570e-21
  3 KSP Residual norm 2.986369469883e-09
Residual norms for fieldsplit_FE_split_ solve.
0 KSP Residual norm 7.813223137340e-16
1 KSP Residual norm 4.793103235095e-18
2 KSP Residual norm 1.615526128222e-18
3 KSP Residual norm 1.179102504397e-18
4 KSP Residual norm 1.427467627551e-19
5 KSP Residual norm 5.177440470993e-20
6 KSP Residual norm 4.703763659148e-20
7 KSP Residual norm 

[petsc-users] Using PCFIELDSPLIT with -pc_fieldsplit_type schur

[David Knezevic]

I have a definite block 2x2 system and I figured it'd be good to apply the
PCFIELDSPLIT functionality with Schur complement, as described in Section
4.5 of the manual.

The A00 block of my matrix is very small so I figured I'd specify a direct
solver (i.e. MUMPS) for that block.

So I did the following:
- PCFieldSplitSetIS to specify the indices of the two splits
- PCFieldSplitGetSubKSP to get the two KSP objects, and to set the solver
and PC types for each (MUMPS for A00, ILU+CG for A11)
- I set -pc_fieldsplit_schur_fact_type full

Below I have pasted the output of "-ksp_view -ksp_monitor -log_view" for a
test case. It seems to converge well, but I'm concerned about the speed
(about 90 seconds, vs. about 1 second if I use a direct solver for the
entire system). I just wanted to check if I'm setting this up in a good way?

Many thanks,
David

---

  0 KSP Residual norm 5.405774214400e+04
  1 KSP Residual norm 1.849649014371e+02
  2 KSP Residual norm 7.462775074989e-02
  3 KSP Residual norm 2.680497175260e-04
KSP Object: 1 MPI processes
  type: cg
  maximum iterations=1000
  tolerances:  relative=1e-06, absolute=1e-50, divergence=1.
  left preconditioning
  using nonzero initial guess
  using PRECONDITIONED norm type for convergence test
PC Object: 1 MPI processes
  type: fieldsplit
FieldSplit with Schur preconditioner, factorization FULL
Preconditioner for the Schur complement formed from A11
Split info:
Split number 0 Defined by IS
Split number 1 Defined by IS
KSP solver for A00 block
  KSP Object:  (fieldsplit_RB_split_)   1 MPI processes
type: preonly
maximum iterations=1, initial guess is zero
tolerances:  relative=1e-05, absolute=1e-50, divergence=1.
left preconditioning
using NONE norm type for convergence test
  PC Object:  (fieldsplit_RB_split_)   1 MPI processes
type: cholesky
  Cholesky: out-of-place factorization
  tolerance for zero pivot 2.22045e-14
  matrix ordering: natural
  factor fill ratio given 0., needed 0.
Factored matrix follows:
  Mat Object:   1 MPI processes
type: seqaij
rows=324, cols=324
package used to perform factorization: mumps
total: nonzeros=3042, allocated nonzeros=3042
total number of mallocs used during MatSetValues calls =0
  MUMPS run parameters:
SYM (matrix type):   2
PAR (host participation):1
ICNTL(1) (output for error): 6
ICNTL(2) (output of diagnostic msg): 0
ICNTL(3) (output for global info):   0
ICNTL(4) (level of printing):0
ICNTL(5) (input mat struct): 0
ICNTL(6) (matrix prescaling):7
ICNTL(7) (sequentia matrix ordering):7
ICNTL(8) (scalling strategy):77
ICNTL(10) (max num of refinements):  0
ICNTL(11) (error analysis):  0
ICNTL(12) (efficiency control):
0
ICNTL(13) (efficiency control):
0
ICNTL(14) (percentage of estimated workspace increase):
20
ICNTL(18) (input mat struct):
0
ICNTL(19) (Shur complement info):
0
ICNTL(20) (rhs sparse pattern):
0
ICNTL(21) (solution struct):
 0
ICNTL(22) (in-core/out-of-core facility):
0
ICNTL(23) (max size of memory can be allocated
locally):0
ICNTL(24) (detection of null pivot rows):
0
ICNTL(25) (computation of a null space basis):
 0
ICNTL(26) (Schur options for rhs or solution):
 0
ICNTL(27) (experimental parameter):
-24
ICNTL(28) (use parallel or sequential ordering):
 1
ICNTL(29) (parallel ordering):
 0
ICNTL(30) (user-specified set of entries in inv(A)):
 0
ICNTL(31) (factors is discarded in the solve phase):
 0
ICNTL(33) (compute determinant):
 0
CNTL(1) (relative pivoting threshold):  0.01
CNTL(2) (stopping criterion of refinement): 1.49012e-08
CNTL(3) (absolute pivoting threshold):  0.
CNTL(4) (value of static pivoting): -1.
CNTL(5) (fixation for null pivots): 0.
RINFO(1) (local estimated flops for the elimination
after analysis):
  [0] 29394.
RINFO(2) (local estimated flops for the assembly after
factorization):
 

Re: [petsc-users] Leveraging IIoT & Moving to Zero Unplanned Downtime

[David Knezevic]

Apologies, the email below was not meant to be sent to this list.

David


On Wed, Jan 4, 2017 at 7:39 AM, Thomas Leurent  wrote:

> Akselos News
> View this email in your browser
> 
>
> Dear All,
>
> I invite you to read my new post
> 
>  on
> how coupling the Industrial Internet of Things (IIoT) with physics-based
> digital twins will be the path to materialise immense savings from the
> industry current IIoT investments.
>
> Upcoming touch points: if you're interested in discussing more this month, 
> I'll
> be attending WEF Annual Meeting
> 
> in Davos, January 17-20 and MIT-ILP
> 
>  event
> in Japan, January 27th. If you are attending, I hope to see you there.
>
> Thomas,
>
> 
>
> If you work in the Oil & Gas industry, you’ve probably noticed that its
> future is being shaped around the notion of Zero Unplanned Downtime. The
> offshore sector, in particular, is plagued with unplanned downtime that is
> costing it billions in unseen value. And the problem is only going to get
> worse: as operating expenses are scaled down, capital projects are cut or
> deferred, and aging assets are operated beyond their original design life. 
> Read
> more of this post.
> 
>
>
> 
> Twitter
> 
>
> 
> LinkedIn
> 
>  Email 
>
> 
> YouTube
> 
>
> *Copyright ©  2016 Akselos, All rights reserved.*
> You are on this list because you are affiliated to Akselos
>
> *Our mailing address is:*
> Akselos, EPFL Innovation Park, Building D
> Lausanne 1015, Switzerland
> Add us to your address book
> 
>  update your preferences
> 
> or unsubscribe from this list
> 
>
> [image: Email Marketing Powered by MailChimp]
> 
>

Re: [petsc-users] Issue updating MUMPS ictnl after failed solve

[David Knezevic]

On Tue, Sep 20, 2016 at 2:24 PM, Hong  wrote:

> David:
> I'll patch petsc-maint (v3.7), then merge it to petsc-master.
> It might take 1-2 days for regression tests.
>
>
I pulled the patched version and it works well for me now, thanks!

David




> On Tue, Sep 20, 2016 at 2:01 PM, Hong  wrote:
>>
>>> David :
>>> This is a bug in PETSc. A change
>>> $ git diff ../../../pc/impls/factor/cholesky/cholesky.c
>>> diff --git a/src/ksp/pc/impls/factor/cholesky/cholesky.c
>>> b/src/ksp/pc/impls/factor/cholesky/cholesky.c
>>> index 953d551..cc28369 100644
>>> --- a/src/ksp/pc/impls/factor/cholesky/cholesky.c
>>> +++ b/src/ksp/pc/impls/factor/cholesky/cholesky.c
>>> @@ -141,9 +141,7 @@ static PetscErrorCode PCSetUp_Cholesky(PC pc)
>>>
>>>  ierr = MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&(
>>> (PC_Factor*)dir)->info);CHKERRQ(ierr);
>>>  ierr = MatFactorGetError(((PC_Factor*)dir)->fact,);CHKERRQ(ierr
>>> );
>>> -if (err) { /* FactorNumeric() fails */
>>> -  pc->failedreason = (PCFailedReason)err;
>>> -}
>>> +pc->failedreason = (PCFailedReason)err;
>>>}
>>>
>>> fixed the problem. I'll fix this problem in petsc-release, including
>>> other routines.
>>> Thanks for reporting the bug and sending matrix.dat. Let us know
>>> whenever you encounter problem using PETSc.
>>>
>>
>>
>> OK, great, thanks for the fix.
>>
>> Will this fix be included in the next patch release of 3.7?
>>
>> Thanks,
>> David
>>
>>
>>
>>
>
>

Re: [petsc-users] Issue updating MUMPS ictnl after failed solve

[David Knezevic]

On Tue, Sep 20, 2016 at 2:01 PM, Hong  wrote:

> David :
> This is a bug in PETSc. A change
> $ git diff ../../../pc/impls/factor/cholesky/cholesky.c
> diff --git a/src/ksp/pc/impls/factor/cholesky/cholesky.c
> b/src/ksp/pc/impls/factor/cholesky/cholesky.c
> index 953d551..cc28369 100644
> --- a/src/ksp/pc/impls/factor/cholesky/cholesky.c
> +++ b/src/ksp/pc/impls/factor/cholesky/cholesky.c
> @@ -141,9 +141,7 @@ static PetscErrorCode PCSetUp_Cholesky(PC pc)
>
>  ierr = MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&(
> (PC_Factor*)dir)->info);CHKERRQ(ierr);
>  ierr = MatFactorGetError(((PC_Factor*)dir)->fact,);CHKERRQ(ierr);
> -if (err) { /* FactorNumeric() fails */
> -  pc->failedreason = (PCFailedReason)err;
> -}
> +pc->failedreason = (PCFailedReason)err;
>}
>
> fixed the problem. I'll fix this problem in petsc-release, including other
> routines.
> Thanks for reporting the bug and sending matrix.dat. Let us know whenever
> you encounter problem using PETSc.
>


OK, great, thanks for the fix.

Will this fix be included in the next patch release of 3.7?

Thanks,
David

Re: [petsc-users] Issue updating MUMPS ictnl after failed solve

[David Knezevic]

On Mon, Sep 19, 2016 at 11:04 PM, Hong  wrote:

> David :
> I did following:
>
> PC  pc;
> Mat F;
> ierr = KSPGetPC(ksp,);CHKERRQ(ierr);
> ierr = PCReset(pc);CHKERRQ(ierr);
> ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
> ierr = PCSetType(pc,PCCHOLESKY);CHKERRQ(ierr);
>
> ierr = PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);CHKERRQ(
> ierr);
> ierr = PCFactorSetUpMatSolverPackage(pc);CHKERRQ(ierr);
> ierr = PCFactorGetMatrix(pc,);CHKERRQ(ierr);
> ierr = MatMumpsSetIcntl(F,14,30);CHKERRQ(ierr);
>
> ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
>
> Then it resolves the matrix equation with ICNTL(14)=30.
> Attached is modified petsc/src/ksp/ksp/examples/tutorials/ex10.c.
> Using in with your matrix.dat, I get
>
> mpiexec -n 4 ./ex10 -f0 matrix.dat -rhs 0 -ksp_reason
> Number of iterations =   0
> KSPConvergedReason: -11
>  Reset PC with ICNTL(14)=30 ...
> KSPConvergedReason: 2
>

Hong,

Thanks very much for your test code. I get the same output as you when I
run "mpiexec -n 4 ./ex10 -f0 matrix.dat -rhs 0 -ksp_reason".

However, I used KSPPREONLY in my original test code, and if I
add KSPSetType(ksp,KSPPREONLY) in your modified exc10.c after the line
KSPCreate(PETSC_COMM_WORLD,) then I get the following output:

mpiexec -np 4 ./mumps_test-opt -f0 matrix.dat -rhs 0 -ksp_reason
Number of iterations =   0
KSPConvergedReason: -11
 Reset PC with ICNTL(14)=30 ...
KSPConvergedReason: -11

So it seems like the icntl data is not being updated when we use PREONLY.
Do you know how to fix this?

Thanks,
David

Re: [petsc-users] Issue updating MUMPS ictnl after failed solve

[David Knezevic]

On Mon, Sep 19, 2016 at 9:45 PM, Fande Kong <fdkong...@gmail.com> wrote:

> Placing PCReset(PC pc) before the second kspsolve might works.
>
> Fande Kong,
>
> On Mon, Sep 19, 2016 at 7:38 PM, murat keçeli <kec...@gmail.com> wrote:
>
>> Another guess: maybe you also need KSPSetUp(ksp); before the second
>> KSPSolve(ksp,b,x);.
>>
>> Murat Keceli
>>
>
Thanks for the suggestions. I just tried these, and they didn't work either
unfortunately.

David





> ​
>>
>> On Mon, Sep 19, 2016 at 8:33 PM, David Knezevic <
>> david.kneze...@akselos.com> wrote:
>>
>>> On Mon, Sep 19, 2016 at 7:26 PM, Dave May <dave.mayhe...@gmail.com>
>>> wrote:
>>>
>>>>
>>>>
>>>> On 19 September 2016 at 21:05, David Knezevic <
>>>> david.kneze...@akselos.com> wrote:
>>>>
>>>>> When I use MUMPS via PETSc, one issue is that it can sometimes fail
>>>>> with MUMPS error -9, which means that MUMPS didn't allocate a big enough
>>>>> workspace. This can typically be fixed by increasing MUMPS icntl 14, e.g.
>>>>> via the command line option -mat_mumps_icntl_14.
>>>>>
>>>>> However, instead of having to run several times with different command
>>>>> line options, I'd like to be able to automatically increment icntl 14 
>>>>> value
>>>>> in a loop until the solve succeeds.
>>>>>
>>>>> I have a saved matrix which fails when I use it for a solve with MUMPS
>>>>> with 4 MPI processes and the default ictnl values, so I'm using this to
>>>>> check that I can achieve the automatic icntl 14 update, as described 
>>>>> above.
>>>>> (The matrix is 14MB so I haven't attached it here, but I'd be happy to 
>>>>> send
>>>>> it to anyone else who wants to try this test case out.)
>>>>>
>>>>> I've pasted some test code below which provides a simple test of this
>>>>> idea using two solves. The first solve uses the default value of icntl 14,
>>>>> which fails, and then we update icntl 14 to 30 and solve again. The second
>>>>> solve should succeed since icntl 14 of 30 is sufficient for MUMPS to
>>>>> succeed in this case, but for some reason the second solve still fails.
>>>>>
>>>>> Below I've also pasted the output from -ksp_view, and you can see that
>>>>> ictnl 14 is being updated correctly (see the ICNTL(14) lines in the
>>>>> output), so it's not clear to me why the second solve fails. It seems like
>>>>> MUMPS is ignoring the update to the ictnl value?
>>>>>
>>>>
>>>> I believe this parameter is utilized during the numerical factorization
>>>> phase.
>>>> In your code, the operator hasn't changed, however you haven't
>>>> signalled to the KSP that you want to re-perform the numerical
>>>> factorization.
>>>> You can do this by calling KSPSetOperators() before your second solve.
>>>> I think if you do this (please try it), the factorization will be
>>>> performed again and the new value of icntl will have an effect.
>>>>
>>>> Note this is a wild stab in the dark - I haven't dug through the
>>>> petsc-mumps code in detail...
>>>>
>>>
>>> That sounds like a plausible guess to me, but unfortunately it didn't
>>> work. I added KSPSetOperators(ksp,A,A); before the second solve and I
>>> got the same behavior as before.
>>>
>>> Thanks,
>>> David
>>>
>>>
>>>
>>>
>>>
>>>> 
>>>>> -
>>>>> Test code:
>>>>>
>>>>>   Mat A;
>>>>>   MatCreate(PETSC_COMM_WORLD,);
>>>>>   MatSetType(A,MATMPIAIJ);
>>>>>
>>>>>   PetscViewer petsc_viewer;
>>>>>   PetscViewerBinaryOpen( PETSC_COMM_WORLD,
>>>>>  "matrix.dat",
>>>>>  FILE_MODE_READ,
>>>>>  _viewer);
>>>>>   MatLoad(A, petsc_viewer);
>>>>>   PetscViewerDestroy(_viewer);
>>>>>
>>>>>   PetscInt m, n;
>>>>>   MatGetSize(A, , );
>>>>>
>>>>>   Vec x;
>>>>>   VecCreate(PET

Re: [petsc-users] Issue updating MUMPS ictnl after failed solve

[David Knezevic]

On Mon, Sep 19, 2016 at 7:26 PM, Dave May <dave.mayhe...@gmail.com> wrote:

>
>
> On 19 September 2016 at 21:05, David Knezevic <david.kneze...@akselos.com>
> wrote:
>
>> When I use MUMPS via PETSc, one issue is that it can sometimes fail with
>> MUMPS error -9, which means that MUMPS didn't allocate a big enough
>> workspace. This can typically be fixed by increasing MUMPS icntl 14, e.g.
>> via the command line option -mat_mumps_icntl_14.
>>
>> However, instead of having to run several times with different command
>> line options, I'd like to be able to automatically increment icntl 14 value
>> in a loop until the solve succeeds.
>>
>> I have a saved matrix which fails when I use it for a solve with MUMPS
>> with 4 MPI processes and the default ictnl values, so I'm using this to
>> check that I can achieve the automatic icntl 14 update, as described above.
>> (The matrix is 14MB so I haven't attached it here, but I'd be happy to send
>> it to anyone else who wants to try this test case out.)
>>
>> I've pasted some test code below which provides a simple test of this
>> idea using two solves. The first solve uses the default value of icntl 14,
>> which fails, and then we update icntl 14 to 30 and solve again. The second
>> solve should succeed since icntl 14 of 30 is sufficient for MUMPS to
>> succeed in this case, but for some reason the second solve still fails.
>>
>> Below I've also pasted the output from -ksp_view, and you can see that
>> ictnl 14 is being updated correctly (see the ICNTL(14) lines in the
>> output), so it's not clear to me why the second solve fails. It seems like
>> MUMPS is ignoring the update to the ictnl value?
>>
>
> I believe this parameter is utilized during the numerical factorization
> phase.
> In your code, the operator hasn't changed, however you haven't signalled
> to the KSP that you want to re-perform the numerical factorization.
> You can do this by calling KSPSetOperators() before your second solve.
> I think if you do this (please try it), the factorization will be
> performed again and the new value of icntl will have an effect.
>
> Note this is a wild stab in the dark - I haven't dug through the
> petsc-mumps code in detail...
>

That sounds like a plausible guess to me, but unfortunately it didn't work.
I added KSPSetOperators(ksp,A,A); before the second solve and I got the
same behavior as before.

Thanks,
David





> 
>> -
>> Test code:
>>
>>   Mat A;
>>   MatCreate(PETSC_COMM_WORLD,);
>>   MatSetType(A,MATMPIAIJ);
>>
>>   PetscViewer petsc_viewer;
>>   PetscViewerBinaryOpen( PETSC_COMM_WORLD,
>>  "matrix.dat",
>>  FILE_MODE_READ,
>>  _viewer);
>>   MatLoad(A, petsc_viewer);
>>   PetscViewerDestroy(_viewer);
>>
>>   PetscInt m, n;
>>   MatGetSize(A, , );
>>
>>   Vec x;
>>   VecCreate(PETSC_COMM_WORLD,);
>>   VecSetSizes(x,PETSC_DECIDE,m);
>>   VecSetFromOptions(x);
>>   VecSet(x,1.0);
>>
>>   Vec b;
>>   VecDuplicate(x,);
>>
>>   KSP ksp;
>>   PC pc;
>>
>>   KSPCreate(PETSC_COMM_WORLD,);
>>   KSPSetOperators(ksp,A,A);
>>
>>   KSPSetType(ksp,KSPPREONLY);
>>   KSPGetPC(ksp,);
>>
>>   PCSetType(pc,PCCHOLESKY);
>>
>>   PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
>>   PCFactorSetUpMatSolverPackage(pc);
>>
>>   KSPSetFromOptions(ksp);
>>   KSPSetUp(ksp);
>>
>>   KSPSolve(ksp,b,x);
>>
>>   {
>> KSPConvergedReason reason;
>> KSPGetConvergedReason(ksp, );
>> std::cout << "converged reason: " << reason << std::endl;
>>   }
>>
>>   Mat F;
>>   PCFactorGetMatrix(pc,);
>>   MatMumpsSetIcntl(F,14,30);
>>
>>   KSPSolve(ksp,b,x);
>>
>>   {
>> KSPConvergedReason reason;
>> KSPGetConvergedReason(ksp, );
>> std::cout << "converged reason: " << reason << std::endl;
>>   }
>>
>> 
>> -
>> -ksp_view output (ICNTL(14) changes from 20 to 30, but we get "converged
>> reason: -11" for both solves)
>>
>> KSP Object: 4 MPI processes
>>   type: preonly
>>   maximum iterations=1, initial guess is zero
>>   tolerances:  relative=1e-05, absolute=1e-50, divergence=1.
&g

[petsc-users] Issue updating MUMPS ictnl after failed solve

[David Knezevic]

When I use MUMPS via PETSc, one issue is that it can sometimes fail with
MUMPS error -9, which means that MUMPS didn't allocate a big enough
workspace. This can typically be fixed by increasing MUMPS icntl 14, e.g.
via the command line option -mat_mumps_icntl_14.

However, instead of having to run several times with different command line
options, I'd like to be able to automatically increment icntl 14 value in a
loop until the solve succeeds.

I have a saved matrix which fails when I use it for a solve with MUMPS with
4 MPI processes and the default ictnl values, so I'm using this to check
that I can achieve the automatic icntl 14 update, as described above. (The
matrix is 14MB so I haven't attached it here, but I'd be happy to send it
to anyone else who wants to try this test case out.)

I've pasted some test code below which provides a simple test of this idea
using two solves. The first solve uses the default value of icntl 14, which
fails, and then we update icntl 14 to 30 and solve again. The second solve
should succeed since icntl 14 of 30 is sufficient for MUMPS to succeed in
this case, but for some reason the second solve still fails.

Below I've also pasted the output from -ksp_view, and you can see that
ictnl 14 is being updated correctly (see the ICNTL(14) lines in the
output), so it's not clear to me why the second solve fails. It seems like
MUMPS is ignoring the update to the ictnl value?

Thanks,
David


-
Test code:

  Mat A;
  MatCreate(PETSC_COMM_WORLD,);
  MatSetType(A,MATMPIAIJ);

  PetscViewer petsc_viewer;
  PetscViewerBinaryOpen( PETSC_COMM_WORLD,
 "matrix.dat",
 FILE_MODE_READ,
 _viewer);
  MatLoad(A, petsc_viewer);
  PetscViewerDestroy(_viewer);

  PetscInt m, n;
  MatGetSize(A, , );

  Vec x;
  VecCreate(PETSC_COMM_WORLD,);
  VecSetSizes(x,PETSC_DECIDE,m);
  VecSetFromOptions(x);
  VecSet(x,1.0);

  Vec b;
  VecDuplicate(x,);

  KSP ksp;
  PC pc;

  KSPCreate(PETSC_COMM_WORLD,);
  KSPSetOperators(ksp,A,A);

  KSPSetType(ksp,KSPPREONLY);
  KSPGetPC(ksp,);

  PCSetType(pc,PCCHOLESKY);

  PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
  PCFactorSetUpMatSolverPackage(pc);

  KSPSetFromOptions(ksp);
  KSPSetUp(ksp);

  KSPSolve(ksp,b,x);

  {
KSPConvergedReason reason;
KSPGetConvergedReason(ksp, );
std::cout << "converged reason: " << reason << std::endl;
  }

  Mat F;
  PCFactorGetMatrix(pc,);
  MatMumpsSetIcntl(F,14,30);

  KSPSolve(ksp,b,x);

  {
KSPConvergedReason reason;
KSPGetConvergedReason(ksp, );
std::cout << "converged reason: " << reason << std::endl;
  }


-
-ksp_view output (ICNTL(14) changes from 20 to 30, but we get "converged
reason: -11" for both solves)

KSP Object: 4 MPI processes
  type: preonly
  maximum iterations=1, initial guess is zero
  tolerances:  relative=1e-05, absolute=1e-50, divergence=1.
  left preconditioning
  using NONE norm type for convergence test
PC Object: 4 MPI processes
  type: cholesky
Cholesky: out-of-place factorization
tolerance for zero pivot 2.22045e-14
matrix ordering: natural
factor fill ratio given 0., needed 0.
  Factored matrix follows:
Mat Object: 4 MPI processes
  type: mpiaij
  rows=22878, cols=22878
  package used to perform factorization: mumps
  total: nonzeros=3361617, allocated nonzeros=3361617
  total number of mallocs used during MatSetValues calls =0
MUMPS run parameters:
  SYM (matrix type):   2
  PAR (host participation):1
  ICNTL(1) (output for error): 6
  ICNTL(2) (output of diagnostic msg): 0
  ICNTL(3) (output for global info):   0
  ICNTL(4) (level of printing):0
  ICNTL(5) (input mat struct): 0
  ICNTL(6) (matrix prescaling):7
  ICNTL(7) (sequentia matrix ordering):7
  ICNTL(8) (scalling strategy):77
  ICNTL(10) (max num of refinements):  0
  ICNTL(11) (error analysis):  0
  ICNTL(12) (efficiency control): 0
  ICNTL(13) (efficiency control): 0
  ICNTL(14) (percentage of estimated workspace increase): 20
  ICNTL(18) (input mat struct):   3
  ICNTL(19) (Shur complement info):   0
  ICNTL(20) (rhs sparse pattern): 0
  ICNTL(21) (solution struct):1
  ICNTL(22) (in-core/out-of-core facility):   0
  ICNTL(23) (max size of memory can be allocated locally):0
   

Re: [petsc-users] MUMPS error reporting in PETSc-3.7

[David Knezevic]

On Mon, Jun 6, 2016 at 4:16 PM, Barry Smith <bsm...@mcs.anl.gov> wrote:

>
>   David,
>
>Actually all the routines for accessing the MUMPS information already
> exist. I have added this to the manual page for MATSOLVERMUMPS
>
> Notes: When a MUMPS factorization fails inside a KSP solve, for
> example with a KSP_DIVERGED_PCSETUP_FAILED, one can find the MUMPS
> information about the failure by calling
> $  KSPGetPC(ksp,);
> $  PCFactorGetMatrix(pc,);
> $  MatMumpsGetInfo(mat,);
> $  MatMumpsGetInfog(mat,); etc.
>
> .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage,
> MatMumpsSetICntl(), MatMumpsGetIcntl(), MatMumpsSetCntl(),
> MatMumpsGetCntl(), MatMumpsGetInfo(), MatMumpsGetInfog(),
> MatMumpsGetRinfo(), MatMumpsGetRinfog(), KSPGetPC(), PCGetFactor(),
> PCFactorGetMatrix()
>
> Please let us know if this does not work
>


Thanks! That looks good. I've changed my code to use this, and I'll let you
know if I have any issues with it.

Thanks,
David




> On Jun 5, 2016, at 12:45 PM, David Knezevic <david.kneze...@akselos.com>
> wrote:
> >
> > On Sun, Jun 5, 2016 at 1:42 PM, Barry Smith <bsm...@mcs.anl.gov> wrote:
> >
> > > On Jun 5, 2016, at 12:12 PM, David Knezevic <
> david.kneze...@akselos.com> wrote:
> > >
> > > In PETSc-3.6, when I hit an error using the MUMPS solver, I would get
> an error like this:
> > >
> > > [0]PETSC ERROR: - Error Message
> --
> > > [0]PETSC ERROR: Error in external library
> > > [0]PETSC ERROR: Error reported by MUMPS in numerical factorization
> phase: INFO(1)=-9, INFO(2)=351513
> > >
> > > In PETSc-3.7, the solver now reports "KSP_DIVERGED_PCSETUP_FAILED"
> instead. The new behavior is definitely better, so thanks for the update!
> However, I'd still be curious to be able to see what error codes MUMPS
> returned (which was previously provided in INFO(1) and INFO(2)), since that
> can be helpful for diagnosing the reason for the solver failure. Is that
> info still available somehow?
> >
> >David,
> >
> >  There are a couple of ways to obtain them both clunky, we should
> also provide a simple function call interface such as
> MatMumpsGetInfo(Mat,int *infog,int *info);
> >
> >  Anyways right now you can run with -info and it prints the
> information (along with too much other information) or you can run with
> -ksp_error_if_not_converged and then it reverts back to the old behavior of
> stopping immediately and printing the error message and info values.
> >
> >
> > OK, got it, thanks!
> >
> > David
>
>

Re: [petsc-users] MUMPS error reporting in PETSc-3.7

[David Knezevic]

On Sun, Jun 5, 2016 at 1:42 PM, Barry Smith <bsm...@mcs.anl.gov> wrote:

>
> > On Jun 5, 2016, at 12:12 PM, David Knezevic <david.kneze...@akselos.com>
> wrote:
> >
> > In PETSc-3.6, when I hit an error using the MUMPS solver, I would get an
> error like this:
> >
> > [0]PETSC ERROR: - Error Message
> --
> > [0]PETSC ERROR: Error in external library
> > [0]PETSC ERROR: Error reported by MUMPS in numerical factorization
> phase: INFO(1)=-9, INFO(2)=351513
> >
> > In PETSc-3.7, the solver now reports "KSP_DIVERGED_PCSETUP_FAILED"
> instead. The new behavior is definitely better, so thanks for the update!
> However, I'd still be curious to be able to see what error codes MUMPS
> returned (which was previously provided in INFO(1) and INFO(2)), since that
> can be helpful for diagnosing the reason for the solver failure. Is that
> info still available somehow?
>
>David,
>
>  There are a couple of ways to obtain them both clunky, we should also
> provide a simple function call interface such as MatMumpsGetInfo(Mat,int
> *infog,int *info);
>
>  Anyways right now you can run with -info and it prints the
> information (along with too much other information) or you can run with
> -ksp_error_if_not_converged and then it reverts back to the old behavior of
> stopping immediately and printing the error message and info values.
>


OK, got it, thanks!

David

[petsc-users] MUMPS error reporting in PETSc-3.7

[David Knezevic]

In PETSc-3.6, when I hit an error using the MUMPS solver, I would get an
error like this:

[0]PETSC ERROR: - Error Message
--
[0]PETSC ERROR: Error in external library
[0]PETSC ERROR: Error reported by MUMPS in numerical factorization phase:
INFO(1)=-9, INFO(2)=351513

In PETSc-3.7, the solver now reports "KSP_DIVERGED_PCSETUP_FAILED" instead.
The new behavior is definitely better, so thanks for the update! However,
I'd still be curious to be able to see what error codes MUMPS returned
(which was previously provided in INFO(1) and INFO(2)), since that can be
helpful for diagnosing the reason for the solver failure. Is that info
still available somehow?

Thanks,
David

Re: [petsc-users] eigensolution error with slepc

[David Knezevic]

On Mon, Apr 4, 2016 at 1:27 PM, Jose E. Roman <jro...@dsic.upv.es> wrote:

>
> > El 4 abr 2016, a las 19:13, David Knezevic <david.kneze...@akselos.com>
> escribió:
> >
> > OK, thanks, I'll have a look at that paper.
> >
> > And just to confirm that I've understood properly: You're saying that
> the GNHEP solver should converge more robustly than the GHIEP solver for
> symmetric indefinite problems? Are there any advantages to using GHIEP over
> GNHEP (e.g. is it faster?), or do you just recommend to always use GNHEP
> for symmetric indefinite problems?
> >
> > Apologies if the paper already answers those questions, I'll read
> through it soon.
> >
> > Thanks,
> > David
>
> In terms of linear algebra, a symmetric-indefinite matrix pencil is not
> symmetric (no guarantee it has real eigenvalues and a full set of linearly
> independent eigenvectors). So in principle a non-symmetric algorithm should
> be used. The GHIEP solver tries to enforce (pseudo-) symmetry throughout
> the computation. This may provide more accurate results or maybe converge
> faster, but the cost per iteration is higher (since it uses B-inner
> products rather than standard inner products as in Arnoldi). Maybe in
> problems arising in buckling analysis of structures it is always safe to
> use GHIEP, I don't know.



OK, got it, thanks!

David

Re: [petsc-users] eigensolution error with slepc

[David Knezevic]

On Mon, Apr 4, 2016 at 12:47 PM, Jose E. Roman <jro...@dsic.upv.es> wrote:

>
> > El 4 abr 2016, a las 18:32, David Knezevic <david.kneze...@akselos.com>
> escribió:
> >
> > Hi Jose,
> >
> > I'm interested to know more about this comment from your email below:
> > "You can also try the symmetric-indefinite solver, but this is not
> recommended since it is not numerically stable."
> >
> > I frequently use SLEPc's symmetric-indefinite solver for buckling
> eigenvalue problems (which are naturally symmetric and indefinite), and
> this works well in my experience. But, based on your comment, I was
> wondering if you would recommend a different solver for this case?
> >
> > Thanks,
> > David
>
> The symmetric-indefinite solver is based on non-symmetric Lanczos, which
> is not guaranteed to be numerically stable (it does not rely on orthogonal
> projections as Arnoldi). In practice, it means that the method may break
> down or get unstable (in our implementation we detect the latter). For many
> problems the indefinite solver will provide the solution without problems,
> but no guarantees it will always work.
>
> Details of the indefinite solver can be found in section 3 of this paper:
> http://dx.doi.org/10.1007/s10543-016-0601-5
>
> Jose
>
>
OK, thanks, I'll have a look at that paper.

And just to confirm that I've understood properly: You're saying that the
GNHEP solver should converge more robustly than the GHIEP solver for
symmetric indefinite problems? Are there any advantages to using GHIEP over
GNHEP (e.g. is it faster?), or do you just recommend to always use GNHEP
for symmetric indefinite problems?

Apologies if the paper already answers those questions, I'll read through
it soon.

Thanks,
David

Re: [petsc-users] eigensolution error with slepc

[David Knezevic]

Hi Jose,

I'm interested to know more about this comment from your email below:
"You can also try the symmetric-indefinite solver, but this is not
recommended since it is not numerically stable."

I frequently use SLEPc's symmetric-indefinite solver for buckling
eigenvalue problems (which are naturally symmetric and indefinite), and
this works well in my experience. But, based on your comment, I was
wondering if you would recommend a different solver for this case?

Thanks,
David



On Mon, Apr 4, 2016 at 12:18 PM, Jose E. Roman  wrote:

>
> > El 4 abr 2016, a las 18:06, Manav Bhatia 
> escribió:
> >
> > Hi Jose,
> >
> >I also read these matrices into matlab and found the eigenvalues as
> >
> > >>A = PetscBinaryRead('A.petsc’);
> > >>B = PetscBinaryRead(‘B.petsc’);
> > >> [v,d] = eigs(A,B)
> > (*** got a lot of output about poor-conditioning ***)
> > >> diag(d)
> >
> > ans =
> >
> >1.0e-05 *
> >
> >-0.2219
> > 0.0229
> > 0.0229
> > 0.0025
> > 0.0019
> > 0.0014
> >
> > >> So, one of these is turning out to be negative, but I am still
> getting numbers.
> >
> >
> >
> > Using these matrices with ex7 produces an error:
> >
> > Dhcp-90-243:T_400_V_0 manav$
> /Users/manav/Documents/codes/numerical_lib/slepc/slepc-3.6.2/src/eps/examples/tutorials/ex7
> -f1 A.petsc -f2 B.petsc -eps_gen_hermitian -eps_view -eps_nev 1
> >
> > Generalized eigenproblem stored in file.
> >
> >  Reading REAL matrices from binary files...
> > [0]PETSC ERROR: - Error Message
> --
> > [0]PETSC ERROR: Error in external library
> > [0]PETSC ERROR: Error in Lapack xSTEQR 15
> > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
> for trouble shooting.
> > [0]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015
> > [0]PETSC ERROR:
> /Users/manav/Documents/codes/numerical_lib/slepc/slepc-3.6.2/src/eps/examples/tutorials/ex7
> on a arch-darwin-cxx-opt named ws243-49.walker.dynamic.msstate.edu by
> manav Mon Apr  4 11:05:30 2016
> > [0]PETSC ERROR: Configure options
> --prefix=/Users/manav/Documents/codes/numerical_lib/petsc/petsc-3.6.3/../
> --CC=mpicc --CXX=mpicxx --FC=mpif90 --with-clanguage=c++ --with-fortran=0
> --with-mpiexec=/opt/local/bin/mpiexec --with-shared-libraries=1 --with-x=1
> --with-x-dir=/opt/X11 --with-debugging=0
> --with-lapack-lib=/usr/lib/liblapack.dylib
> --with-blas-lib=/usr/lib/libblas.dylib --download-superlu=yes
> --download-superlu_dist=yes --download-suitesparse=yes --download-mumps=yes
> --download-scalapack=yes --download-parmetis=yes
> --download-parmetis-shared=1 --download-metis=yes --download-metis-shared=1
> --download-hypre=yes --download-hypre-shared=1 --download-ml=yes
> --download-ml-shared=1 --download-sundials=yes --download-sundials-shared=1
> > [0]PETSC ERROR: #1 DSSolve_HEP_QR() line 495 in
> /Users/manav/Documents/codes/numerical_lib/slepc/slepc-3.6.2/src/sys/classes/ds/impls/hep/dshep.c
> > [0]PETSC ERROR: #2 DSSolve() line 543 in
> /Users/manav/Documents/codes/numerical_lib/slepc/slepc-3.6.2/src/sys/classes/ds/interface/dsops.c
> > [0]PETSC ERROR: #3 EPSSolve_KrylovSchur_Symm() line 68 in
> /Users/manav/Documents/codes/numerical_lib/slepc/slepc-3.6.2/src/eps/impls/krylov/krylovschur/ks-symm.c
> > [0]PETSC ERROR: #4 EPSSolve() line 101 in
> /Users/manav/Documents/codes/numerical_lib/slepc/slepc-3.6.2/src/eps/interface/epssolve.c
> > [0]PETSC ERROR: #5 main() line 147 in
> /Users/manav/Documents/codes/numerical_lib/slepc/slepc-3.6.2/src/eps/examples/tutorials/ex7.c
> > [0]PETSC ERROR: PETSc Option Table entries:
> > [0]PETSC ERROR: -eps_gen_hermitian
> > [0]PETSC ERROR: -eps_nev 1
> > [0]PETSC ERROR: -eps_view
> > [0]PETSC ERROR: -f1 A.petsc
> > [0]PETSC ERROR: -f2 B.petsc
> > [0]PETSC ERROR: -matload_block_size 1
> > [0]PETSC ERROR: End of Error Message ---send entire
> error message to petsc-ma...@mcs.anl.gov--
> >
> --
> > MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
> > with errorcode 76.
> >
> > NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
> > You may or may not see output from other processes, depending on
> > exactly when Open MPI kills them.
> >
> --
> >
>
> The problem is that now K is indefinite. The check for indefinite B-matrix
> was not done correctly. This is fixed in 3.6.3, where you should get a more
> informative error:
>
> [0]PETSC ERROR: The inner product is not well defined: indefinite matrix
>
> When considered as a non-symmetric problem, the solver returns reasonable
> output:
>
> $ ./ex7 -f1 ~/tmp/bhatia/A.petsc -f2 ~/tmp/bhatia/B.petsc
>
> Generalized eigenproblem stored in file.
>
>  Reading REAL matrices from binary files...
>  Number of iterations of the method: 1
>  Number of linear 

Re: [petsc-users] PCFactorSetUpMatSolverPackage with SNES

[David Knezevic]

Hi Barry,

I have added an error checker so the code will produce a useful error
> message instead of just crashing here in the code.
>

This sounds helpful. I assume this is in the dev branch? I'm using 3.6.1,
but I gather from this list that 3.7 will be out soon, so I'll switch to
that once it's available.



> You can only call this routine AFTER the matrix has been provided to
> the KSP object. Normally with SNES this won't happen until after the solve
> has started (hence is not easy for you to put in your parameters) but
> adding the code I suggest below might work
>


OK, makes sense. Unfortunately the change below didn't help, I still get
the same segfault. But it's not a big deal, since I went with Matt's
suggestion which seems to be working well.

Thanks,
David




> > On Mar 1, 2016, at 12:52 PM, David Knezevic <david.kneze...@akselos.com>
> wrote:
> >
> > Based on KSP ex52, I use PCFactorSetUpMatSolverPackage in the process of
> setting various MUMPS ictnl options. This works fine for me when I'm
> solving linear problems.
> >
> > I then wanted to use PCFactorSetUpMatSolverPackage with the PC from a
> SNES object. I tried to do this with the following code (after calling
> SNESCreate, SNESSetFunction, and SNESSetJacobian):
> >
> > KSP snes_ksp;
> > SNESGetKSP(snes, _ksp);
> KSPSetOperators(ksp,mat,pmat);
>  /* mat and pmat (which may be the same) are the Jacobian arguments to
> SNESSetJacobian(), the matrices must exist and types set but they don't
> have to have the correct Jacobian entries at this point (since the
> nonlinear solver has not started you cannot know the Jacobian entries yet.)
> > PC snes_pc;
> > KSPGetPC(snes_ksp, _pc);
> PCSetType(snes_pc,PCLU);
> > PCFactorSetMatSolverPackage(snes_pc, MATSOLVERMUMPS);
> > PCFactorSetUpMatSolverPackage(snes_pc);
>
> Let me know if this does not work and where it goes wrong.
> >
> > However, I get a segfault on the call to PCFactorSetUpMatSolverPackage
> in this case. I was wondering what I need to do to make this work?
> >
> > Note that I want to set the MUMPS ictnl parameters via code rather than
> via the commandline since sometimes MUMPS fails (e.g. with error -9 due to
> a workspace size that is too small) and I need to automatically re-run the
> solve with different ictnl values when this happens.
> >
> > Thanks,
> > David
> >
> >
>
>

Re: [petsc-users] PCFactorSetUpMatSolverPackage with SNES

[David Knezevic]

On Tue, Mar 1, 2016 at 1:56 PM, Matthew Knepley <knep...@gmail.com> wrote:

> On Tue, Mar 1, 2016 at 12:52 PM, David Knezevic <
> david.kneze...@akselos.com> wrote:
>
>> Based on KSP ex52, I use PCFactorSetUpMatSolverPackage in the process of
>> setting various MUMPS ictnl options. This works fine for me when I'm
>> solving linear problems.
>>
>> I then wanted to use PCFactorSetUpMatSolverPackage with the PC from a
>> SNES object. I tried to do this with the following code (after calling
>> SNESCreate, SNESSetFunction, and SNESSetJacobian):
>>
>> KSP snes_ksp;
>> SNESGetKSP(snes, _ksp);
>> PC snes_pc;
>> KSPGetPC(snes_ksp, _pc);
>> PCFactorSetMatSolverPackage(snes_pc, MATSOLVERMUMPS);
>> PCFactorSetUpMatSolverPackage(snes_pc);
>>
>> However, I get a segfault on the call to PCFactorSetUpMatSolverPackage in
>> this case. I was wondering what I need to do to make this work?
>>
>> Note that I want to set the MUMPS ictnl parameters via code rather than
>> via the commandline since sometimes MUMPS fails (e.g. with error -9 due to
>> a workspace size that is too small) and I need to automatically re-run the
>> solve with different ictnl values when this happens.
>>
>
> That is a good reason. However, I would organize this differently. I would
> still set the type from the command line.
> Then later in your code, after SNES is setup correctly, you get out the PC
> and reset the icntl values if you have a
> failure. Its very difficult to get the setup logic correct, and its one of
> the most error prone parts of PETSc. RIght now,
> the most reliable way to do things is to have all the information
> available up front in options.
>
> Someday, I want to write a small DFA piece in PETSc that can encode all
> this logic so that simple errors people
> make can be diagnosed early with nice error messages.
>
>   Thanks,
>
>  Matt
>

OK, thanks for the info, I'll go with that approach.

David

[petsc-users] PCFactorSetUpMatSolverPackage with SNES

[David Knezevic]

Based on KSP ex52, I use PCFactorSetUpMatSolverPackage in the process of
setting various MUMPS ictnl options. This works fine for me when I'm
solving linear problems.

I then wanted to use PCFactorSetUpMatSolverPackage with the PC from a SNES
object. I tried to do this with the following code (after calling
SNESCreate, SNESSetFunction, and SNESSetJacobian):

KSP snes_ksp;
SNESGetKSP(snes, _ksp);
PC snes_pc;
KSPGetPC(snes_ksp, _pc);
PCFactorSetMatSolverPackage(snes_pc, MATSOLVERMUMPS);
PCFactorSetUpMatSolverPackage(snes_pc);

However, I get a segfault on the call to PCFactorSetUpMatSolverPackage in
this case. I was wondering what I need to do to make this work?

Note that I want to set the MUMPS ictnl parameters via code rather than via
the commandline since sometimes MUMPS fails (e.g. with error -9 due to a
workspace size that is too small) and I need to automatically re-run the
solve with different ictnl values when this happens.

Thanks,
David

Re: [petsc-users] MatPtAP

[David Knezevic]

On Wed, Feb 24, 2016 at 1:06 AM, Barry Smith <bsm...@mcs.anl.gov> wrote:

>
>   Yes, could definitely be faster with a custom code but a bit of a pain
> to write.
>


OK, thanks for your input!

David





>
> > On Feb 23, 2016, at 10:08 PM, David Knezevic <david.kneze...@akselos.com>
> wrote:
> >
> > On Tue, Feb 23, 2016 at 10:40 PM, Barry Smith <bsm...@mcs.anl.gov>
> wrote:
> >
> >   A custom MatPtAP would almost surely pay off, but it is only an
> optimization so you need to ask if this computation is your main "blocker"
> to getting more "science" done.  Can you send a picture of the exact
> structure of P?
> >
> > I've attached a picture of the sparsity pattern of P. There are 0s and
> 1s on the diagonal, and some dense columns in the first few columns.
> >
> > David
> >
> >
> >
> >
> >
> >
> >
> >
> > > On Feb 23, 2016, at 9:35 PM, David Knezevic <
> david.kneze...@akselos.com> wrote:
> > >
> > > I'm using MatPtAP, which works well for me, but in some examples I've
> tested the PtAP calculation dominates the overall solve time (e.g. see
> attached -log_summary output).
> > >
> > > In my case, A is a stiffness matrix, and P is the identity matrix
> except for a small number of columns (e.g. about 10 or so) which are dense.
> > >
> > > In this situation, I was wondering if there is a more efficient way to
> proceed than using MatPtAP? For example, would it be noticeably faster to
> calculate P^T A P directly using MatMults for the dense columns, rather
> than using MatPtAP?
> > >
> > > Thanks!
> > > David
> > >
> > >
> > > 
> >
> >
> > 
>
>

[petsc-users] MatPtAP

[David Knezevic]

I'm using MatPtAP, which works well for me, but in some examples I've
tested the PtAP calculation dominates the overall solve time (e.g. see
attached -log_summary output).

In my case, A is a stiffness matrix, and P is the identity matrix except
for a small number of columns (e.g. about 10 or so) which are dense.

In this situation, I was wondering if there is a more efficient way to
proceed than using MatPtAP? For example, would it be noticeably faster to
calculate P^T A P directly using MatMults for the dense columns, rather
than using MatPtAP?

Thanks!
David
-- PETSc Performance Summary: 
--

/home/dknez/akselos-dev/scrbe/build/bin/fe_solver-opt_real on a 
arch-linux2-c-opt named david-Lenovo with 1 processor, by dknez Tue Feb 23 
14:17:28 2016
Using Petsc Release Version 3.6.1, Jul, 22, 2015 

 Max   Max/MinAvg  Total 
Time (sec):   1.214e+02  1.0   1.214e+02
Objects:  4.150e+02  1.0   4.150e+02
Flops:1.294e+09  1.0   1.294e+09  1.294e+09
Flops/sec:1.066e+07  1.0   1.066e+07  1.066e+07
MPI Messages: 0.000e+00  0.0   0.000e+00  0.000e+00
MPI Message Lengths:  0.000e+00  0.0   0.000e+00  0.000e+00
MPI Reductions:   0.000e+00  0.0

Flop counting convention: 1 flop = 1 real number operation of type 
(multiply/divide/add/subtract)
e.g., VecAXPY() for real vectors of length N --> 2N 
flops
and VecAXPY() for complex vectors of length N --> 
8N flops

Summary of Stages:   - Time --  - Flops -  --- Messages ---  -- 
Message Lengths --  -- Reductions --
Avg %Total Avg %Total   counts   %Total 
Avg %Total   counts   %Total 
 0:  Main Stage: 1.2142e+02 100.0%  1.2945e+09 100.0%  0.000e+00   0.0%  
0.000e+000.0%  0.000e+00   0.0% 


See the 'Profiling' chapter of the users' manual for details on interpreting 
output.
Phase summary info:
   Count: number of times phase was executed
   Time and Flops: Max - maximum over all processors
   Ratio - ratio of maximum to minimum over all processors
   Mess: number of messages sent
   Avg. len: average message length (bytes)
   Reduct: number of global reductions
   Global: entire computation
   Stage: stages of a computation. Set stages with PetscLogStagePush() and 
PetscLogStagePop().
  %T - percent time in this phase %F - percent flops in this phase
  %M - percent messages in this phase %L - percent message lengths in 
this phase
  %R - percent reductions in this phase
   Total Mflop/s: 10e-6 * (sum of flops over all processors)/(max time over all 
processors)

EventCount  Time (sec) Flops
 --- Global ---  --- Stage ---   Total
   Max Ratio  Max Ratio   Max  Ratio  Mess   Avg len Reduct 
 %T %F %M %L %R  %T %F %M %L %R Mflop/s


--- Event Stage 0: Main Stage

VecNorm   15 1.0 2.5439e-04 1.0 1.12e+06 1.0 0.0e+00 0.0e+00 
0.0e+00  0  0  0  0  0   0  0  0  0  0  4384
VecCopy  102 1.0 2.6283e-03 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 
0.0e+00  0  0  0  0  0   0  0  0  0  0 0
VecSet   315 1.0 3.1855e-03 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 
0.0e+00  0  0  0  0  0   0  0  0  0  0 0
VecAXPY  104 1.0 3.3436e-03 1.0 1.25e+06 1.0 0.0e+00 0.0e+00 
0.0e+00  0  0  0  0  0   0  0  0  0  0   374
VecWAXPY   7 1.0 4.2009e-04 1.0 2.60e+05 1.0 0.0e+00 0.0e+00 
0.0e+00  0  0  0  0  0   0  0  0  0  0   620
VecAssemblyBegin 368 1.0 1.2109e-03 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 
0.0e+00  0  0  0  0  0   0  0  0  0  0 0
VecAssemblyEnd   368 1.0 2.5558e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 
0.0e+00  0  0  0  0  0   0  0  0  0  0 0
VecScatterBegin  106 1.0 8.7023e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 
0.0e+00  0  0  0  0  0   0  0  0  0  0 0
VecReduceArith21 1.0 5.3144e-04 1.0 1.56e+06 1.0 0.0e+00 0.0e+00 
0.0e+00  0  0  0  0  0   0  0  0  0  0  2940
VecReduceComm  7 1.0 1.1921e-06 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 
0.0e+00  0  0  0  0  0   0  0  0  0  0 0
MatMultAdd41 1.0 3.2815e-02 1.0 5.99e+07 1.0 0.0e+00 0.0e+00 
0.0e+00  0  5  0  0  0   0  5  0  0  0  1824
MatMultTrAdd   9 1.0 6.2582e-03 1.0 8.70e+06 1.0 0.0e+00 0.0e+00 
0.0e+00  0  1  0  0  0   0  1  0  0  0  1390
MatSolve  19 1.0 7.3646e-01 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 

Re: [petsc-users] Convergence trouble with SNES and CG

[David Knezevic]

On Mon, Feb 22, 2016 at 7:10 PM, Matthew Knepley <knep...@gmail.com> wrote:

> On Mon, Feb 22, 2016 at 6:02 PM, David Knezevic <
> david.kneze...@akselos.com> wrote:
>
>> On Mon, Feb 22, 2016 at 6:57 PM, Barry Smith <bsm...@mcs.anl.gov> wrote:
>>
>>>
>>>   Add -ksp_converged_reason
>>>
>>
>>
>> That gives "DIVERGED_INDEFINITE_PC". The full output is below:
>>
>
> Just for testing, replace CG with GMRES.
>

With GMRES the first four nonlinear iterations behave much that same as
before, converging at KSP iteration 150 or so.

On the 5th nonlinear iteration, GMRES doesn't give a KSP convergence
failure, but the convergence seems to stagnate. It's currently at KSP
iteration 2000, and tolerance of 6e-3 (and still running).

David

Re: [petsc-users] Convergence trouble with SNES and CG

[David Knezevic]

On Mon, Feb 22, 2016 at 6:57 PM, Barry Smith  wrote:

>
>   Add -ksp_converged_reason
>


That gives "DIVERGED_INDEFINITE_PC". The full output is below:


  105 KSP Residual norm 7.099534889984e-05
  106 KSP Residual norm 6.615351528474e-05
  107 KSP Residual norm 6.303646443688e-05
  Linear solve converged due to CONVERGED_RTOL iterations 107
  NL step  5, |residual|_2 = 1.996487e+01
0 KSP Residual norm 5.844456247091e+00
1 KSP Residual norm 1.308756943973e+00
2 KSP Residual norm 7.545133187463e-01
3 KSP Residual norm 1.346068942607e+00
  Linear solve did not converge due to DIVERGED_INDEFINITE_PC iterations 4
SNES Object: 1 MPI processes
  type: newtonls
  maximum iterations=50, maximum function evaluations=1
  tolerances: relative=1e-08, absolute=1e-50, solution=1e-08
  total number of linear solver iterations=535
  total number of function evaluations=11
  norm schedule ALWAYS
  SNESLineSearch Object:   1 MPI processes
type: bt
  interpolation: cubic
  alpha=1.00e-04
maxstep=1.00e+08, minlambda=1.00e-12
tolerances: relative=1.00e-08, absolute=1.00e-15,
lambda=1.00e-08
maximum iterations=40
  KSP Object:   1 MPI processes
type: cg
maximum iterations=1, initial guess is zero
tolerances:  relative=1e-05, absolute=1e-50, divergence=1
left preconditioning
using NATURAL norm type for convergence test
  PC Object:   1 MPI processes
type: ilu
  ILU: out-of-place factorization
  0 levels of fill
  tolerance for zero pivot 2.22045e-14
  matrix ordering: natural
  factor fill ratio given 1, needed 1
Factored matrix follows:
  Mat Object:   1 MPI processes
type: seqaij
rows=967608, cols=967608
package used to perform factorization: petsc
total: nonzeros=7.55539e+07, allocated nonzeros=7.55539e+07
total number of mallocs used during MatSetValues calls =0
  using I-node routines: found 322536 nodes, limit used is 5
linear system matrix = precond matrix:
Mat Object:() 1 MPI processes
  type: seqaij
  rows=967608, cols=967608
  total: nonzeros=7.55539e+07, allocated nonzeros=7.55539e+07
  total number of mallocs used during MatSetValues calls =0
using I-node routines: found 322536 nodes, limit used is 5
Number of nonlinear iterations: 5
Nonlinear solver convergence/divergence reason: DIVERGED_LINEAR_SOLVE

[petsc-users] Convergence trouble with SNES and CG

[David Knezevic]

I'm solving a "linear elasticity with contact" problem using SNES. The
matrix is symmetric, and it corresponds to the standard linear elasticity
matrix, plus some extra nonlinear "penalty spring" terms due to the contact
conditions.

The solve works fine if I use a direct solver (e.g. MUMPS). I've also found
that CG usually works fine too, but in one specific case it fails to
converge. (I thought that penalty springs + iterative solvers may not be a
good combination, but it does seem to work fine in most cases.)

Below I've pasted the (last part of the) output of "-ksp_type cg
-ksp_norm_type natural -snes_view -ksp_monitor" for the case that fails.

If anyone has some suggestions about other solver/preconditioner options
that would be worth trying, that would be most appreciated.

Thanks,
David

---

  102 KSP Residual norm 8.910585826247e-05
  103 KSP Residual norm 8.395178827803e-05
  104 KSP Residual norm 7.758478472168e-05
  105 KSP Residual norm 7.099534889984e-05
  106 KSP Residual norm 6.615351528474e-05
  107 KSP Residual norm 6.303646443688e-05
  NL step  5, |residual|_2 = 1.996487e+01
0 KSP Residual norm 5.844456247091e+00
1 KSP Residual norm 1.308756943973e+00
2 KSP Residual norm 7.545133187463e-01
3 KSP Residual norm 1.346068942607e+00
SNES Object: 1 MPI processes
  type: newtonls
  maximum iterations=50, maximum function evaluations=1
  tolerances: relative=1e-08, absolute=1e-50, solution=1e-08
  total number of linear solver iterations=535
  total number of function evaluations=11
  norm schedule ALWAYS
  SNESLineSearch Object:   1 MPI processes
type: bt
  interpolation: cubic
  alpha=1.00e-04
maxstep=1.00e+08, minlambda=1.00e-12
tolerances: relative=1.00e-08, absolute=1.00e-15,
lambda=1.00e-08
maximum iterations=40
  KSP Object:   1 MPI processes
type: cg
maximum iterations=1, initial guess is zero
tolerances:  relative=1e-05, absolute=1e-50, divergence=1
left preconditioning
using NATURAL norm type for convergence test
  PC Object:   1 MPI processes
type: ilu
  ILU: out-of-place factorization
  0 levels of fill
  tolerance for zero pivot 2.22045e-14
  matrix ordering: natural
  factor fill ratio given 1, needed 1
Factored matrix follows:
  Mat Object:   1 MPI processes
type: seqaij
rows=967608, cols=967608
package used to perform factorization: petsc
total: nonzeros=7.55539e+07, allocated nonzeros=7.55539e+07
total number of mallocs used during MatSetValues calls =0
  using I-node routines: found 322536 nodes, limit used is 5
linear system matrix = precond matrix:
Mat Object:() 1 MPI processes
  type: seqaij
  rows=967608, cols=967608
  total: nonzeros=7.55539e+07, allocated nonzeros=7.55539e+07
  total number of mallocs used during MatSetValues calls =0
using I-node routines: found 322536 nodes, limit used is 5
Number of nonlinear iterations: 5
Nonlinear solver convergence/divergence reason: DIVERGED_LINEAR_SOLVE

[petsc-users] SNES NEWTONLS serial vs. parallel

[David Knezevic]

I'm using NEWTONLS (with mumps for the linear solves) to do a nonlinear PDE
solve. It converges well when I use 1 core. When I use 2 or more cores, the
line search stagnates. I've pasted the output of -snes_linesearch_monitor
below in these two cases.

I was wondering if this implies that I must have a bug in parallel, or if
perhaps the NEWTONLS solver can behave slightly differently in parallel?

Thanks,
David

---

*Serial case:*
  NL step  0, |residual|_2 = 4.714515e-02
  Line search: gnorm after quadratic fit 7.862867755323e-02
  Line search: Cubically determined step, current gnorm
4.663945043239e-02 lambda=1.4276549921126183e-02
  NL step  1, |residual|_2 = 4.663945e-02
  Line search: gnorm after quadratic fit 6.977268575068e-02
  Line search: Cubically determined step, current gnorm
4.594912794004e-02 lambda=2.3644825912085998e-02
  NL step  2, |residual|_2 = 4.594913e-02
  Line search: gnorm after quadratic fit 5.502067932478e-02
  Line search: Cubically determined step, current gnorm
4.494531294405e-02 lambda=4.1260497615261321e-02
  NL step  3, |residual|_2 = 4.494531e-02
  Line search: gnorm after quadratic fit 5.415371063247e-02
  Line search: Cubically determined step, current gnorm
4.392165925471e-02 lambda=3.6375618871780056e-02
  NL step  4, |residual|_2 = 4.392166e-02
  Line search: gnorm after quadratic fit 4.631663976615e-02
  Line search: Cubically determined step, current gnorm
4.246200798775e-02 lambda=5.0003e-02
  NL step  5, |residual|_2 = 4.246201e-02
  Line search: gnorm after quadratic fit 4.222105321728e-02
  Line search: Quadratically determined step,
lambda=1.0001e-01
  NL step  6, |residual|_2 = 4.222105e-02
  Line search: gnorm after quadratic fit 4.026081251872e-02
  Line search: Quadratically determined step,
lambda=1.0001e-01
  NL step  7, |residual|_2 = 4.026081e-02
  Line search: gnorm after quadratic fit 3.776439532346e-02
  Line search: Quadratically determined step,
lambda=1.0001e-01
  NL step  8, |residual|_2 = 3.776440e-02
  Line search: gnorm after quadratic fit 3.659796311121e-02
  Line search: Quadratically determined step,
lambda=1.0001e-01
  NL step  9, |residual|_2 = 3.659796e-02
  Line search: gnorm after quadratic fit 3.423207664901e-02
  Line search: Quadratically determined step,
lambda=1.0001e-01
  NL step 10, |residual|_2 = 3.423208e-02
  Line search: gnorm after quadratic fit 3.116928452225e-02
  Line search: Quadratically determined step,
lambda=1.0001e-01
  NL step 11, |residual|_2 = 3.116928e-02
  Line search: gnorm after quadratic fit 2.874310955274e-02
  Line search: Quadratically determined step,
lambda=1.0001e-01
  NL step 12, |residual|_2 = 2.874311e-02
  Line search: gnorm after quadratic fit 2.587826662305e-02
  Line search: Quadratically determined step,
lambda=1.0001e-01
  NL step 13, |residual|_2 = 2.587827e-02
  Line search: gnorm after quadratic fit 2.344161073075e-02
  Line search: Quadratically determined step,
lambda=1.0001e-01
  NL step 14, |residual|_2 = 2.344161e-02
  Line search: gnorm after quadratic fit 2.187719889554e-02
  Line search: Quadratically determined step,
lambda=1.0001e-01
  NL step 15, |residual|_2 = 2.187720e-02
  Line search: gnorm after quadratic fit 1.983089075086e-02
  Line search: Quadratically determined step,
lambda=1.0001e-01
  NL step 16, |residual|_2 = 1.983089e-02
  Line search: gnorm after quadratic fit 1.791227711151e-02
  Line search: Quadratically determined step,
lambda=1.0001e-01
  NL step 17, |residual|_2 = 1.791228e-02
  Line search: gnorm after quadratic fit 1.613250573900e-02
  Line search: Quadratically determined step,
lambda=1.0001e-01
  NL step 18, |residual|_2 = 1.613251e-02
  Line search: gnorm after quadratic fit 1.455841843183e-02
  Line search: Quadratically determined step,
lambda=1.0001e-01
  NL step 19, |residual|_2 = 1.455842e-02
  Line search: gnorm after quadratic fit 1.321849780208e-02
  Line search: Quadratically determined step,
lambda=1.0574876450981290e-01
  NL step 20, |residual|_2 = 1.321850e-02
  Line search: gnorm after quadratic fit 9.209641609489e-03
  Line search: Quadratically determined step,
lambda=3.0589684959139674e-01
  NL step 21, |residual|_2 = 9.209642e-03
  Line search: gnorm after quadratic fit 7.590942028574e-03
  Line search: Quadratically determined step,
lambda=2.0920305898507460e-01
  NL step 22, |residual|_2 = 7.590942e-03
  Line search: gnorm after quadratic fit 4.373918927227e-03
  Line search: Quadratically determined step,
lambda=4.2379743128074154e-01
  NL step 23, |residual|_2 = 4.373919e-03
  Line 

Re: [petsc-users] SNES NEWTONLS serial vs. parallel

[David Knezevic]

OK, will do, thanks.

David



On Wed, Jan 13, 2016 at 4:05 PM, Barry Smith <bsm...@mcs.anl.gov> wrote:

>
>   Since you are using a direct solver almost for sure a bug in your
> parallel function or parallel Jacobian.
>
>Try -snes_mf_operator   try -snes_fdtry -snes_type test  as three
> different approaches to see what is going on.
>
>Barry
>
> > On Jan 13, 2016, at 2:51 PM, David Knezevic <david.kneze...@akselos.com>
> wrote:
> >
> > Oops! I pasted the wrong text for the serial case. The correct text is
> below:
> >
> > Serial case:
> >   NL step  0, |residual|_2 = 4.714515e-02
> >   Line search: gnorm after quadratic fit 7.862867755130e-02
> >   Line search: Cubically determined step, current gnorm
> 4.663945044088e-02 lambda=1.4276549223307832e-02
> >   NL step  1, |residual|_2 = 4.663945e-02
> >   Line search: gnorm after quadratic fit 6.977268532963e-02
> >   Line search: Cubically determined step, current gnorm
> 4.594912791877e-02 lambda=2.3644826349821228e-02
> >   NL step  2, |residual|_2 = 4.594913e-02
> >   Line search: gnorm after quadratic fit 5.502067915588e-02
> >   Line search: Cubically determined step, current gnorm
> 4.494531287593e-02 lambda=4.1260496881982515e-02
> >   NL step  3, |residual|_2 = 4.494531e-02
> >   Line search: gnorm after quadratic fit 5.415371014813e-02
> >   Line search: Cubically determined step, current gnorm
> 4.392165909219e-02 lambda=3.6375617606865668e-02
> >   NL step  4, |residual|_2 = 4.392166e-02
> >   Line search: gnorm after quadratic fit 4.631663907262e-02
> >   Line search: Cubically determined step, current gnorm
> 4.246200768767e-02 lambda=5.0003e-02
> >   NL step  5, |residual|_2 = 4.246201e-02
> >   Line search: gnorm after quadratic fit 4.222105256158e-02
> >   Line search: Quadratically determined step,
> lambda=1.0001e-01
> >   NL step  6, |residual|_2 = 4.222105e-02
> >   Line search: gnorm after quadratic fit 4.026081168915e-02
> >   Line search: Quadratically determined step,
> lambda=1.0001e-01
> >   NL step  7, |residual|_2 = 4.026081e-02
> >   Line search: gnorm after quadratic fit 3.776439443011e-02
> >   Line search: Quadratically determined step,
> lambda=1.0001e-01
> >   NL step  8, |residual|_2 = 3.776439e-02
> >   Line search: gnorm after quadratic fit 3.659796213553e-02
> >   Line search: Quadratically determined step,
> lambda=1.0001e-01
> >   NL step  9, |residual|_2 = 3.659796e-02
> >   Line search: gnorm after quadratic fit 3.423207563496e-02
> >   Line search: Quadratically determined step,
> lambda=1.0001e-01
> >   NL step 10, |residual|_2 = 3.423208e-02
> >   Line search: gnorm after quadratic fit 3.116928356075e-02
> >   Line search: Quadratically determined step,
> lambda=1.0001e-01
> >   NL step 11, |residual|_2 = 3.116928e-02
> >   Line search: gnorm after quadratic fit 2.874310673331e-02
> >   Line search: Quadratically determined step,
> lambda=1.0001e-01
> >   NL step 12, |residual|_2 = 2.874311e-02
> >   Line search: gnorm after quadratic fit 2.587826447631e-02
> >   Line search: Quadratically determined step,
> lambda=1.0001e-01
> >   NL step 13, |residual|_2 = 2.587826e-02
> >   Line search: gnorm after quadratic fit 2.344160918669e-02
> >   Line search: Quadratically determined step,
> lambda=1.0001e-01
> >   NL step 14, |residual|_2 = 2.344161e-02
> >   Line search: gnorm after quadratic fit 2.187719801063e-02
> >   Line search: Quadratically determined step,
> lambda=1.0001e-01
> >   NL step 15, |residual|_2 = 2.187720e-02
> >   Line search: gnorm after quadratic fit 1.983089025936e-02
> >   Line search: Quadratically determined step,
> lambda=1.0001e-01
> >   NL step 16, |residual|_2 = 1.983089e-02
> >   Line search: gnorm after quadratic fit 1.791227696650e-02
> >   Line search: Quadratically determined step,
> lambda=1.0001e-01
> >   NL step 17, |residual|_2 = 1.791228e-02
> >   Line search: gnorm after quadratic fit 1.613250592206e-02
> >   Line search: Quadratically determined step,
> lambda=1.0001e-01
> >   NL step 18, |residual|_2 = 1.613251e-02
> >   Line search: gnorm after quadratic fit 1.455841890804e-02
> >   Line search: Quadratically determined step,
> lambda=1.0001e-01
> >   NL step 

Re: [petsc-users] SNES NEWTONLS serial vs. parallel

[David Knezevic]

  Line search: Quadratically determined step,
lambda=4.3574150080610474e-01
  NL step 26, |residual|_2 = 1.803192e-03
  Line search: Using full step: fnorm 1.803191839408e-03 gnorm
9.015954497317e-04
  NL step 27, |residual|_2 = 9.015954e-04
  Line search: Using full step: fnorm 9.015954497317e-04 gnorm
1.390181456520e-13
  NL step 28, |residual|_2 = 1.390181e-13
Number of nonlinear iterations: 28


On Wed, Jan 13, 2016 at 3:48 PM, David Knezevic <david.kneze...@akselos.com>
wrote:

> I'm using NEWTONLS (with mumps for the linear solves) to do a nonlinear
> PDE solve. It converges well when I use 1 core. When I use 2 or more cores,
> the line search stagnates. I've pasted the output of
> -snes_linesearch_monitor below in these two cases.
>
> I was wondering if this implies that I must have a bug in parallel, or if
> perhaps the NEWTONLS solver can behave slightly differently in parallel?
>
> Thanks,
> David
>
>
> ---
>
>
>
> *Parallel case:*
>   NL step  0, |residual|_2 = 4.714515e-02
>   Line search: gnorm after quadratic fit 7.862867755323e-02
>   Line search: Cubically determined step, current gnorm
> 4.663945043239e-02 lambda=1.4276549921126183e-02
>   NL step  1, |residual|_2 = 4.663945e-02
>   Line search: gnorm after quadratic fit 6.977268575068e-02
>   Line search: Cubically determined step, current gnorm
> 4.594912794004e-02 lambda=2.3644825912085998e-02
>   NL step  2, |residual|_2 = 4.594913e-02
>   Line search: gnorm after quadratic fit 5.502067932478e-02
>   Line search: Cubically determined step, current gnorm
> 4.494531294405e-02 lambda=4.1260497615261321e-02
>   NL step  3, |residual|_2 = 4.494531e-02
>   Line search: gnorm after quadratic fit 5.415371063247e-02
>   Line search: Cubically determined step, current gnorm
> 4.392165925471e-02 lambda=3.6375618871780056e-02
>   NL step  4, |residual|_2 = 4.392166e-02
>   Line search: gnorm after quadratic fit 4.631663976615e-02
>   Line search: Cubically determined step, current gnorm
> 4.246200798775e-02 lambda=5.0003e-02
>   NL step  5, |residual|_2 = 4.246201e-02
>   Line search: gnorm after quadratic fit 4.222105321728e-02
>   Line search: Quadratically determined step,
> lambda=1.0001e-01
>   NL step  6, |residual|_2 = 4.222105e-02
>   Line search: gnorm after quadratic fit 4.026081251872e-02
>   Line search: Quadratically determined step,
> lambda=1.0001e-01
>   NL step  7, |residual|_2 = 4.026081e-02
>   Line search: gnorm after quadratic fit 3.776439532346e-02
>   Line search: Quadratically determined step,
> lambda=1.0001e-01
>   NL step  8, |residual|_2 = 3.776440e-02
>   Line search: gnorm after quadratic fit 3.659796311121e-02
>   Line search: Quadratically determined step,
> lambda=1.0001e-01
>   NL step  9, |residual|_2 = 3.659796e-02
>   Line search: gnorm after quadratic fit 3.423207664901e-02
>   Line search: Quadratically determined step,
> lambda=1.0001e-01
>   NL step 10, |residual|_2 = 3.423208e-02
>   Line search: gnorm after quadratic fit 3.116928452225e-02
>   Line search: Quadratically determined step,
> lambda=1.0001e-01
>   NL step 11, |residual|_2 = 3.116928e-02
>   Line search: gnorm after quadratic fit 2.874310955274e-02
>   Line search: Quadratically determined step,
> lambda=1.0001e-01
>   NL step 12, |residual|_2 = 2.874311e-02
>   Line search: gnorm after quadratic fit 2.587826662305e-02
>   Line search: Quadratically determined step,
> lambda=1.0001e-01
>   NL step 13, |residual|_2 = 2.587827e-02
>   Line search: gnorm after quadratic fit 2.344161073075e-02
>   Line search: Quadratically determined step,
> lambda=1.0001e-01
>   NL step 14, |residual|_2 = 2.344161e-02
>   Line search: gnorm after quadratic fit 2.187719889554e-02
>   Line search: Quadratically determined step,
> lambda=1.0001e-01
>   NL step 15, |residual|_2 = 2.187720e-02
>   Line search: gnorm after quadratic fit 1.983089075086e-02
>   Line search: Quadratically determined step,
> lambda=1.0001e-01
>   NL step 16, |residual|_2 = 1.983089e-02
>   Line search: gnorm after quadratic fit 1.791227711151e-02
>   Line search: Quadratically determined step,
> lambda=1.0001e-01
>   NL step 17, |residual|_2 = 1.791228e-02
>   Line search: gnorm after quadratic fit 1.613250573900e-02
>   Line search: Quadratically determined step,
> lambda=1.0001e-01
>   NL step 18, |residual|_

Re: [petsc-users] GAMG and zero pivots follow up

[David Knezevic]

To follow up on this, I have run some tests on a smaller test case (same
model, coarser mesh). The options I tried are:

8 MPI processes with:
"-ksp_type gmres -pc_type gamg"
"-ksp_type fgmres -pc_type gamg"
"-ksp_type cg -pc_type gamg"
"-ksp_type cg -pc_type gamg -ksp_norm_type natural"
"-ksp_type cg -pc_type gamg -ksp_norm_type natural -pc_mg_type full"

4 MPI processes with:
"-ksp_type cg -pc_type gamg -ksp_norm_type natural -pc_mg_type full"
"-ksp_type cg -pc_type gamg"
"-ksp_type gmres -pc_type gamg"
"-ksp_type fgmres -pc_type gamg"

Let me know if you'd like me to try any other cases.

The "-ksp_monitor -ksp_view" output is attached. As you can see, all of the
"4 MPI processes" cases worked well, and none of the "8 MPI processes"
cases converged.

Best regards,
David




On Fri, Nov 13, 2015 at 2:48 PM, Matthew Knepley  wrote:

> On Fri, Nov 13, 2015 at 12:28 PM, Jed Brown  wrote:
>
>> Matthew Knepley  writes:
>> > Something very strange is happening here. CG should converge
>> monotonically,
>> > but above it does not. What could be happening?
>>
>> Are you use -ksp_norm_type natural?  CG is not monotone in other norms.
>>
>
> Yikes! I did not check that. Why do we have PRECONDITIONED as the default
> for CG?
>
>Matt
>
> Also, if boundary conditions are enforced using a nonsymmetric
>> formulation (for example), then you can get lack of monotonicity with CG
>> that may not be catastrophic.
>>
>
>
>
> --
> 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
>
--

8 MPI processes with "-ksp_type gmres -pc_type gamg -ksp_max_it 200":

  0 KSP Residual norm 1.243143122667e-03 
  1 KSP Residual norm 6.790564176716e-04 
  2 KSP Residual norm 2.863642167759e-04 
  3 KSP Residual norm 2.169460015654e-04 
  4 KSP Residual norm 1.623823395363e-04 
  5 KSP Residual norm 1.115831456626e-04 
  6 KSP Residual norm 7.905550178072e-05 
  7 KSP Residual norm 6.676506348708e-05 
  8 KSP Residual norm 5.411318965562e-05 
  9 KSP Residual norm 4.946393535811e-05 
 10 KSP Residual norm 4.873373148542e-05 
 11 KSP Residual norm 4.873372878673e-05 
 12 KSP Residual norm 4.832393150132e-05 
 13 KSP Residual norm 4.682644869071e-05 
 14 KSP Residual norm 4.177965048741e-05 
 15 KSP Residual norm 3.561738393315e-05 
 16 KSP Residual norm 3.183178450210e-05 
 17 KSP Residual norm 2.820849441834e-05 
 18 KSP Residual norm 2.411305833934e-05 
 19 KSP Residual norm 2.073531106031e-05 
 20 KSP Residual norm 1.832253875945e-05 
 21 KSP Residual norm 1.613725457732e-05 
 22 KSP Residual norm 1.447115239529e-05 
 23 KSP Residual norm 1.332661204650e-05 
 24 KSP Residual norm 1.248919278483e-05 
 25 KSP Residual norm 1.196188151016e-05 
 26 KSP Residual norm 1.161737695363e-05 
 27 KSP Residual norm 1.153632298642e-05 
 28 KSP Residual norm 1.153553211867e-05 
 29 KSP Residual norm 1.150729646659e-05 
 30 KSP Residual norm 1.134640588584e-05 
 31 KSP Residual norm 1.125651483355e-05 
 32 KSP Residual norm 1.121985823782e-05 
 33 KSP Residual norm 1.121682797994e-05 
 34 KSP Residual norm 1.120526536096e-05 
 35 KSP Residual norm 1.112698441144e-05 
 36 KSP Residual norm 1.099161361515e-05 
 37 KSP Residual norm 1.077160664786e-05 
 38 KSP Residual norm 1.046319447066e-05 
 39 KSP Residual norm 1.002732866634e-05 
 40 KSP Residual norm 9.687406818053e-06 
 41 KSP Residual norm 9.291736292845e-06 
 42 KSP Residual norm 8.787280517217e-06 
 43 KSP Residual norm 8.323595238657e-06 
 44 KSP Residual norm 7.891080867185e-06 
 45 KSP Residual norm 7.537064831605e-06 
 46 KSP Residual norm 7.316511381129e-06 
 47 KSP Residual norm 7.185951262668e-06 
 48 KSP Residual norm 7.117216131634e-06 
 49 KSP Residual norm 7.104770082988e-06 
 50 KSP Residual norm 7.099139305978e-06 
 51 KSP Residual norm 7.038487040610e-06 
 52 KSP Residual norm 6.861458611935e-06 
 53 KSP Residual norm 6.625293689513e-06 
 54 KSP Residual norm 6.329429663991e-06 
 55 KSP Residual norm 5.919056214664e-06 
 56 KSP Residual norm 5.507182558921e-06 
 57 KSP Residual norm 5.196803884679e-06 
 58 KSP Residual norm 5.002188092285e-06 
 59 KSP Residual norm 4.880759791404e-06 
 60 KSP Residual norm 4.770150595855e-06 
 61 KSP Residual norm 4.724469615685e-06 
 62 KSP Residual norm 4.673829760077e-06 
 63 KSP Residual norm 4.629705280910e-06 
 64 KSP Residual norm 4.601474765626e-06 
 65 KSP Residual norm 4.593132269745e-06 
 66 KSP Residual norm 4.593013889961e-06 
 67 KSP Residual norm 4.587628601477e-06 
 68 KSP Residual norm 4.552820908762e-06 
 69 KSP Residual norm 4.477855982146e-06 
 70 KSP Residual norm 4.405304333703e-06 
 71 KSP Residual norm 4.330447444642e-06 
 72 KSP Residual norm 4.237398563528e-06 
 73 KSP Residual norm 4.138174613148e-06 
 74 KSP Residual norm 

Re: [petsc-users] GAMG and zero pivots follow up

[David Knezevic]

On Fri, Nov 13, 2015 at 9:45 AM, Matthew Knepley <knep...@gmail.com> wrote:

> On Thu, Nov 12, 2015 at 10:02 PM, David Knezevic <
> david.kneze...@akselos.com> wrote:
>
>> On Thu, Nov 12, 2015 at 9:58 AM, David Knezevic <
>> david.kneze...@akselos.com> wrote:
>>
>>> On Thu, Nov 12, 2015 at 9:50 AM, Mark Adams <mfad...@lbl.gov> wrote:
>>>
>>>> Note, I suspect the zero pivot is coming from a coarse grid.  I don't
>>>> know why no-inode fixed it.  N.B., this might not be deterministic.
>>>>
>>>> If you run with -info and grep on 'GAMG' you will see the block sizes.
>>>> They should be 3 on the fine grid and the 6 on the coarse grids if you set
>>>> everything up correctly.
>>>>
>>>
>>> OK, thanks for the info, I'll look into this further.
>>>
>>> Though note that I got it to work now without needing no-inode now. The
>>> change I made was to make sure that I matched the order of function calls
>>> from the PETSc GAMG examples. The libMesh code I was using was doing things
>>> somewhat out of order, apparently.
>>>
>>
>>
>> As I mentioned above, GAMG seems to be working fine for me now. Thanks
>> for the help on this.
>>
>> However, I wanted to ask about another 3D elasticity test case I ran in
>> which GAMG didn't converge. The "-ksp_view -ksp_monitor" output is below.
>> Any insights into what options I should set in order to get it to converge
>> in this case would be appreciated. (By the way, ML worked in this case with
>> the default options.)
>>
>> Thanks,
>> David
>>
>> 
>>
>>   0 KSP Residual norm 2.475892877233e-03
>>   1 KSP Residual norm 6.181785698923e-05
>>   2 KSP Residual norm 9.439050821656e-04
>>
>
> Something very strange is happening here. CG should converge
> monotonically, but above it does not. What could be happening?
>
>   a) We could have lost orthogonality
>
>   This seems really unlikely in 2 iterates.
>
>   b) The preconditioner could be nonlinear
>
>  In order to test this, could you run with both GMRES and FGMRES? If
> the convergence is different, then something is wrong
>  with the preconditioner.
>
>  I have looked below, and it seems like this should be linear. Mark,
> could something be happening in GAMG?
>
> Also, you should not be using CG/V-cycle. You should be using 1 iterate of
> FMG, -pc_mg_type full.  I have just been going over
> this in class. It is normal for people to solve way past discretization
> error, but that does not make much sense. Its not hard to
> get a sense of discretization error using a manufactured solution.
>


OK, thanks for this info. I'll do some tests based on your comments above
when I get some time and I'll let you know what happens.

Thanks,
David




> KSP Object: 8 MPI processes
>>   type: cg
>>   maximum iterations=5000
>>   tolerances:  relative=1e-12, absolute=1e-50, divergence=1
>>   left preconditioning
>>   using nonzero initial guess
>>   using PRECONDITIONED norm type for convergence test
>> PC Object: 8 MPI processes
>>   type: gamg
>> MG: type is MULTIPLICATIVE, levels=5 cycles=v
>>   Cycles per PCApply=1
>>   Using Galerkin computed coarse grid matrices
>>   GAMG specific options
>> Threshold for dropping small values from graph 0
>> AGG specific options
>>   Symmetric graph false
>>   Coarse grid solver -- level ---
>> KSP Object:(mg_coarse_) 8 MPI processes
>>   type: gmres
>> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt
>> Orthogonalization with no iterative refinement
>> GMRES: happy breakdown tolerance 1e-30
>>   maximum iterations=1, initial guess is zero
>>   tolerances:  relative=1e-05, absolute=1e-50, divergence=1
>>   left preconditioning
>>   using NONE norm type for convergence test
>> PC Object:(mg_coarse_) 8 MPI processes
>>   type: bjacobi
>> block Jacobi: number of blocks = 8
>> Local solve is same for all blocks, in the following KSP and PC
>> objects:
>>   KSP Object:  (mg_coarse_sub_)   1 MPI processes
>> type: preonly
>> maximum iterations=1, initial guess is zero
>> tolerances:  relative=1e-05, absolute=1e-50, divergence=1
>> left preconditioning
>> using NONE norm type for conver

Re: [petsc-users] GAMG and zero pivots follow up

[David Knezevic]

On Thu, Nov 12, 2015 at 9:58 AM, David Knezevic <david.kneze...@akselos.com>
wrote:

> On Thu, Nov 12, 2015 at 9:50 AM, Mark Adams <mfad...@lbl.gov> wrote:
>
>> Note, I suspect the zero pivot is coming from a coarse grid.  I don't
>> know why no-inode fixed it.  N.B., this might not be deterministic.
>>
>> If you run with -info and grep on 'GAMG' you will see the block sizes.
>> They should be 3 on the fine grid and the 6 on the coarse grids if you set
>> everything up correctly.
>>
>
> OK, thanks for the info, I'll look into this further.
>
> Though note that I got it to work now without needing no-inode now. The
> change I made was to make sure that I matched the order of function calls
> from the PETSc GAMG examples. The libMesh code I was using was doing things
> somewhat out of order, apparently.
>


As I mentioned above, GAMG seems to be working fine for me now. Thanks for
the help on this.

However, I wanted to ask about another 3D elasticity test case I ran in
which GAMG didn't converge. The "-ksp_view -ksp_monitor" output is below.
Any insights into what options I should set in order to get it to converge
in this case would be appreciated. (By the way, ML worked in this case with
the default options.)

Thanks,
David



  0 KSP Residual norm 2.475892877233e-03
  1 KSP Residual norm 6.181785698923e-05
  2 KSP Residual norm 9.439050821656e-04
KSP Object: 8 MPI processes
  type: cg
  maximum iterations=5000
  tolerances:  relative=1e-12, absolute=1e-50, divergence=1
  left preconditioning
  using nonzero initial guess
  using PRECONDITIONED norm type for convergence test
PC Object: 8 MPI processes
  type: gamg
MG: type is MULTIPLICATIVE, levels=5 cycles=v
  Cycles per PCApply=1
  Using Galerkin computed coarse grid matrices
  GAMG specific options
Threshold for dropping small values from graph 0
AGG specific options
  Symmetric graph false
  Coarse grid solver -- level ---
KSP Object:(mg_coarse_) 8 MPI processes
  type: gmres
GMRES: restart=30, using Classical (unmodified) Gram-Schmidt
Orthogonalization with no iterative refinement
GMRES: happy breakdown tolerance 1e-30
  maximum iterations=1, initial guess is zero
  tolerances:  relative=1e-05, absolute=1e-50, divergence=1
  left preconditioning
  using NONE norm type for convergence test
PC Object:(mg_coarse_) 8 MPI processes
  type: bjacobi
block Jacobi: number of blocks = 8
Local solve is same for all blocks, in the following KSP and PC
objects:
  KSP Object:  (mg_coarse_sub_)   1 MPI processes
type: preonly
maximum iterations=1, initial guess is zero
tolerances:  relative=1e-05, absolute=1e-50, divergence=1
left preconditioning
using NONE norm type for convergence test
  PC Object:  (mg_coarse_sub_)   1 MPI processes
type: lu
  LU: out-of-place factorization
  tolerance for zero pivot 2.22045e-14
  using diagonal shift on blocks to prevent zero pivot [INBLOCKS]
  matrix ordering: nd
  factor fill ratio given 5, needed 1
Factored matrix follows:
  Mat Object:   1 MPI processes
type: seqaij
rows=6, cols=6, bs=6
package used to perform factorization: petsc
total: nonzeros=36, allocated nonzeros=36
total number of mallocs used during MatSetValues calls =0
  using I-node routines: found 2 nodes, limit used is 5
linear system matrix = precond matrix:
Mat Object: 1 MPI processes
  type: seqaij
  rows=6, cols=6, bs=6
  total: nonzeros=36, allocated nonzeros=36
  total number of mallocs used during MatSetValues calls =0
using I-node routines: found 2 nodes, limit used is 5
  linear system matrix = precond matrix:
  Mat Object:   8 MPI processes
type: mpiaij
rows=6, cols=6, bs=6
total: nonzeros=36, allocated nonzeros=36
total number of mallocs used during MatSetValues calls =0
  using I-node (on process 0) routines: found 2 nodes, limit used
is 5
  Down solver (pre-smoother) on level 1 ---
KSP Object:(mg_levels_1_) 8 MPI processes
  type: chebyshev
Chebyshev: eigenvalue estimates:  min = 0.146922, max = 1.61615
Chebyshev: eigenvalues estimated using gmres with translations  [0
0.1; 0 1.1]
KSP Object:(mg_levels_1_esteig_) 8 MPI processes
  type: gmres
GMRES: restart=30, using Classical (unmodified) Gram-Schmidt
Orthogonalization with no iterative refinement
GMRES: happy breakdown tolerance 1

Re: [petsc-users] GAMG and zero pivots follow up

[David Knezevic]

On Thu, Nov 12, 2015 at 9:50 AM, Mark Adams <mfad...@lbl.gov> wrote:

> Note, I suspect the zero pivot is coming from a coarse grid.  I don't know
> why no-inode fixed it.  N.B., this might not be deterministic.
>
> If you run with -info and grep on 'GAMG' you will see the block sizes.
> They should be 3 on the fine grid and the 6 on the coarse grids if you set
> everything up correctly.
>

OK, thanks for the info, I'll look into this further.

Though note that I got it to work now without needing no-inode now. The
change I made was to make sure that I matched the order of function calls
from the PETSc GAMG examples. The libMesh code I was using was doing things
somewhat out of order, apparently.

David





On Wed, Nov 11, 2015 at 1:36 PM, David Knezevic <david.kneze...@akselos.com>
wrote:

> On Wed, Nov 11, 2015 at 12:57 PM, David Knezevic <
> david.kneze...@akselos.com> wrote:
>
>> On Wed, Nov 11, 2015 at 12:24 PM, David Knezevic <
>> david.kneze...@akselos.com> wrote:
>>
>>> On Tue, Nov 10, 2015 at 10:28 PM, David Knezevic <
>>> david.kneze...@akselos.com> wrote:
>>>
>>>> On Tue, Nov 10, 2015 at 10:24 PM, Matthew Knepley <knep...@gmail.com>
>>>> wrote:
>>>>
>>>>> On Tue, Nov 10, 2015 at 9:21 PM, David Knezevic <
>>>>> david.kneze...@akselos.com> wrote:
>>>>>
>>>>>> On Tue, Nov 10, 2015 at 10:00 PM, Matthew Knepley <knep...@gmail.com>
>>>>>> wrote:
>>>>>>
>>>>>>> On Tue, Nov 10, 2015 at 8:39 PM, David Knezevic <
>>>>>>> david.kneze...@akselos.com> wrote:
>>>>>>>
>>>>>>>> I'm looking into using GAMG, so I wanted to start with a simple 3D
>>>>>>>> elasticity problem. When I first tried this, I got the following "zero
>>>>>>>> pivot" error:
>>>>>>>>
>>>>>>>>
>>>>>>>> ---
>>>>>>>>
>>>>>>>> [0]PETSC ERROR: Zero pivot in LU factorization:
>>>>>>>> http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot
>>>>>>>> [0]PETSC ERROR: Zero pivot, row 3
>>>>>>>> [0]PETSC ERROR: See
>>>>>>>> http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble
>>>>>>>> shooting.
>>>>>>>> [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015
>>>>>>>> [0]PETSC ERROR:
>>>>>>>> /home/dknez/akselos-dev/scrbe/build/bin/fe_solver-opt_real on a
>>>>>>>> arch-linux2-c-opt named david-Lenovo by dknez Tue Nov 10 21:26:39 2015
>>>>>>>> [0]PETSC ERROR: Configure options --with-shared-libraries=1
>>>>>>>> --with-debugging=0 --download-suitesparse --download-parmetis
>>>>>>>> --download-blacs
>>>>>>>> --with-blas-lapack-dir=/opt/intel/system_studio_2015.2.050/mkl
>>>>>>>> --CXXFLAGS=-Wl,--no-as-needed --download-scalapack --download-mumps
>>>>>>>> --download-metis --download-superlu_dist
>>>>>>>> --prefix=/home/dknez/software/libmesh_install/opt_real/petsc
>>>>>>>> --download-hypre --download-ml
>>>>>>>> [0]PETSC ERROR: #1 PetscKernel_A_gets_inverse_A_5() line 48 in
>>>>>>>> /home/dknez/software/petsc-3.6.1/src/mat/impls/baij/seq/dgefa5.c
>>>>>>>> [0]PETSC ERROR: #2 MatSOR_SeqAIJ_Inode() line 2808 in
>>>>>>>> /home/dknez/software/petsc-3.6.1/src/mat/impls/aij/seq/inode.c
>>>>>>>> [0]PETSC ERROR: #3 MatSOR() line 3697 in
>>>>>>>> /home/dknez/software/petsc-3.6.1/src/mat/interface/matrix.c
>>>>>>>> [0]PETSC ERROR: #4 PCApply_SOR() line 37 in
>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/sor/sor.c
>>>>>>>> [0]PETSC ERROR: #5 PCApply() line 482 in
>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c
>>>>>>>> [0]PETSC ERROR: #6 KSP_PCApply() line 242 in
>>>>>>>> /home/dknez/software/petsc-3.6.1/include/petsc/private/kspimpl.h
>>>>>>>> [0]PETSC ERROR: #7 KSPInitialResidual() line 63 in
>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itres.c
>>>>>>

Re: [petsc-users] GAMG and zero pivots follow up

[David Knezevic]

On Wed, Nov 11, 2015 at 12:24 PM, David Knezevic <david.kneze...@akselos.com
> wrote:

> On Tue, Nov 10, 2015 at 10:28 PM, David Knezevic <
> david.kneze...@akselos.com> wrote:
>
>> On Tue, Nov 10, 2015 at 10:24 PM, Matthew Knepley <knep...@gmail.com>
>> wrote:
>>
>>> On Tue, Nov 10, 2015 at 9:21 PM, David Knezevic <
>>> david.kneze...@akselos.com> wrote:
>>>
>>>> On Tue, Nov 10, 2015 at 10:00 PM, Matthew Knepley <knep...@gmail.com>
>>>> wrote:
>>>>
>>>>> On Tue, Nov 10, 2015 at 8:39 PM, David Knezevic <
>>>>> david.kneze...@akselos.com> wrote:
>>>>>
>>>>>> I'm looking into using GAMG, so I wanted to start with a simple 3D
>>>>>> elasticity problem. When I first tried this, I got the following "zero
>>>>>> pivot" error:
>>>>>>
>>>>>>
>>>>>> ---
>>>>>>
>>>>>> [0]PETSC ERROR: Zero pivot in LU factorization:
>>>>>> http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot
>>>>>> [0]PETSC ERROR: Zero pivot, row 3
>>>>>> [0]PETSC ERROR: See
>>>>>> http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble
>>>>>> shooting.
>>>>>> [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015
>>>>>> [0]PETSC ERROR:
>>>>>> /home/dknez/akselos-dev/scrbe/build/bin/fe_solver-opt_real on a
>>>>>> arch-linux2-c-opt named david-Lenovo by dknez Tue Nov 10 21:26:39 2015
>>>>>> [0]PETSC ERROR: Configure options --with-shared-libraries=1
>>>>>> --with-debugging=0 --download-suitesparse --download-parmetis
>>>>>> --download-blacs
>>>>>> --with-blas-lapack-dir=/opt/intel/system_studio_2015.2.050/mkl
>>>>>> --CXXFLAGS=-Wl,--no-as-needed --download-scalapack --download-mumps
>>>>>> --download-metis --download-superlu_dist
>>>>>> --prefix=/home/dknez/software/libmesh_install/opt_real/petsc
>>>>>> --download-hypre --download-ml
>>>>>> [0]PETSC ERROR: #1 PetscKernel_A_gets_inverse_A_5() line 48 in
>>>>>> /home/dknez/software/petsc-3.6.1/src/mat/impls/baij/seq/dgefa5.c
>>>>>> [0]PETSC ERROR: #2 MatSOR_SeqAIJ_Inode() line 2808 in
>>>>>> /home/dknez/software/petsc-3.6.1/src/mat/impls/aij/seq/inode.c
>>>>>> [0]PETSC ERROR: #3 MatSOR() line 3697 in
>>>>>> /home/dknez/software/petsc-3.6.1/src/mat/interface/matrix.c
>>>>>> [0]PETSC ERROR: #4 PCApply_SOR() line 37 in
>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/sor/sor.c
>>>>>> [0]PETSC ERROR: #5 PCApply() line 482 in
>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c
>>>>>> [0]PETSC ERROR: #6 KSP_PCApply() line 242 in
>>>>>> /home/dknez/software/petsc-3.6.1/include/petsc/private/kspimpl.h
>>>>>> [0]PETSC ERROR: #7 KSPInitialResidual() line 63 in
>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itres.c
>>>>>> [0]PETSC ERROR: #8 KSPSolve_GMRES() line 235 in
>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/gmres/gmres.c
>>>>>> [0]PETSC ERROR: #9 KSPSolve() line 604 in
>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
>>>>>> [0]PETSC ERROR: #10 KSPSolve_Chebyshev() line 381 in
>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/cheby/cheby.c
>>>>>> [0]PETSC ERROR: #11 KSPSolve() line 604 in
>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
>>>>>> [0]PETSC ERROR: #12 PCMGMCycle_Private() line 19 in
>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
>>>>>> [0]PETSC ERROR: #13 PCMGMCycle_Private() line 48 in
>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
>>>>>> [0]PETSC ERROR: #14 PCApply_MG() line 338 in
>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
>>>>>> [0]PETSC ERROR: #15 PCApply() line 482 in
>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c
>>>>>> [0]PETSC ERROR: #16 KSP_PCApply() line 242 in
>>>>>> /home/dknez/software/petsc-3.6.1/include/petsc/private/kspimp

Re: [petsc-users] GAMG and zero pivots follow up

[David Knezevic]

On Tue, Nov 10, 2015 at 10:28 PM, David Knezevic <david.kneze...@akselos.com
> wrote:

> On Tue, Nov 10, 2015 at 10:24 PM, Matthew Knepley <knep...@gmail.com>
> wrote:
>
>> On Tue, Nov 10, 2015 at 9:21 PM, David Knezevic <
>> david.kneze...@akselos.com> wrote:
>>
>>> On Tue, Nov 10, 2015 at 10:00 PM, Matthew Knepley <knep...@gmail.com>
>>> wrote:
>>>
>>>> On Tue, Nov 10, 2015 at 8:39 PM, David Knezevic <
>>>> david.kneze...@akselos.com> wrote:
>>>>
>>>>> I'm looking into using GAMG, so I wanted to start with a simple 3D
>>>>> elasticity problem. When I first tried this, I got the following "zero
>>>>> pivot" error:
>>>>>
>>>>> ---
>>>>>
>>>>> [0]PETSC ERROR: Zero pivot in LU factorization:
>>>>> http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot
>>>>> [0]PETSC ERROR: Zero pivot, row 3
>>>>> [0]PETSC ERROR: See
>>>>> http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble
>>>>> shooting.
>>>>> [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015
>>>>> [0]PETSC ERROR:
>>>>> /home/dknez/akselos-dev/scrbe/build/bin/fe_solver-opt_real on a
>>>>> arch-linux2-c-opt named david-Lenovo by dknez Tue Nov 10 21:26:39 2015
>>>>> [0]PETSC ERROR: Configure options --with-shared-libraries=1
>>>>> --with-debugging=0 --download-suitesparse --download-parmetis
>>>>> --download-blacs
>>>>> --with-blas-lapack-dir=/opt/intel/system_studio_2015.2.050/mkl
>>>>> --CXXFLAGS=-Wl,--no-as-needed --download-scalapack --download-mumps
>>>>> --download-metis --download-superlu_dist
>>>>> --prefix=/home/dknez/software/libmesh_install/opt_real/petsc
>>>>> --download-hypre --download-ml
>>>>> [0]PETSC ERROR: #1 PetscKernel_A_gets_inverse_A_5() line 48 in
>>>>> /home/dknez/software/petsc-3.6.1/src/mat/impls/baij/seq/dgefa5.c
>>>>> [0]PETSC ERROR: #2 MatSOR_SeqAIJ_Inode() line 2808 in
>>>>> /home/dknez/software/petsc-3.6.1/src/mat/impls/aij/seq/inode.c
>>>>> [0]PETSC ERROR: #3 MatSOR() line 3697 in
>>>>> /home/dknez/software/petsc-3.6.1/src/mat/interface/matrix.c
>>>>> [0]PETSC ERROR: #4 PCApply_SOR() line 37 in
>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/sor/sor.c
>>>>> [0]PETSC ERROR: #5 PCApply() line 482 in
>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c
>>>>> [0]PETSC ERROR: #6 KSP_PCApply() line 242 in
>>>>> /home/dknez/software/petsc-3.6.1/include/petsc/private/kspimpl.h
>>>>> [0]PETSC ERROR: #7 KSPInitialResidual() line 63 in
>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itres.c
>>>>> [0]PETSC ERROR: #8 KSPSolve_GMRES() line 235 in
>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/gmres/gmres.c
>>>>> [0]PETSC ERROR: #9 KSPSolve() line 604 in
>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
>>>>> [0]PETSC ERROR: #10 KSPSolve_Chebyshev() line 381 in
>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/cheby/cheby.c
>>>>> [0]PETSC ERROR: #11 KSPSolve() line 604 in
>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
>>>>> [0]PETSC ERROR: #12 PCMGMCycle_Private() line 19 in
>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
>>>>> [0]PETSC ERROR: #13 PCMGMCycle_Private() line 48 in
>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
>>>>> [0]PETSC ERROR: #14 PCApply_MG() line 338 in
>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
>>>>> [0]PETSC ERROR: #15 PCApply() line 482 in
>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c
>>>>> [0]PETSC ERROR: #16 KSP_PCApply() line 242 in
>>>>> /home/dknez/software/petsc-3.6.1/include/petsc/private/kspimpl.h
>>>>> [0]PETSC ERROR: #17 KSPSolve_CG() line 139 in
>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/cg/cg.c
>>>>> [0]PETSC ERROR: #18 KSPSolve() line 604 in
>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
>>>>>
>>>>> 

Re: [petsc-users] GAMG and zero pivots follow up

[David Knezevic]

On Wed, Nov 11, 2015 at 12:57 PM, David Knezevic <david.kneze...@akselos.com
> wrote:

> On Wed, Nov 11, 2015 at 12:24 PM, David Knezevic <
> david.kneze...@akselos.com> wrote:
>
>> On Tue, Nov 10, 2015 at 10:28 PM, David Knezevic <
>> david.kneze...@akselos.com> wrote:
>>
>>> On Tue, Nov 10, 2015 at 10:24 PM, Matthew Knepley <knep...@gmail.com>
>>> wrote:
>>>
>>>> On Tue, Nov 10, 2015 at 9:21 PM, David Knezevic <
>>>> david.kneze...@akselos.com> wrote:
>>>>
>>>>> On Tue, Nov 10, 2015 at 10:00 PM, Matthew Knepley <knep...@gmail.com>
>>>>> wrote:
>>>>>
>>>>>> On Tue, Nov 10, 2015 at 8:39 PM, David Knezevic <
>>>>>> david.kneze...@akselos.com> wrote:
>>>>>>
>>>>>>> I'm looking into using GAMG, so I wanted to start with a simple 3D
>>>>>>> elasticity problem. When I first tried this, I got the following "zero
>>>>>>> pivot" error:
>>>>>>>
>>>>>>>
>>>>>>> ---
>>>>>>>
>>>>>>> [0]PETSC ERROR: Zero pivot in LU factorization:
>>>>>>> http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot
>>>>>>> [0]PETSC ERROR: Zero pivot, row 3
>>>>>>> [0]PETSC ERROR: See
>>>>>>> http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble
>>>>>>> shooting.
>>>>>>> [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015
>>>>>>> [0]PETSC ERROR:
>>>>>>> /home/dknez/akselos-dev/scrbe/build/bin/fe_solver-opt_real on a
>>>>>>> arch-linux2-c-opt named david-Lenovo by dknez Tue Nov 10 21:26:39 2015
>>>>>>> [0]PETSC ERROR: Configure options --with-shared-libraries=1
>>>>>>> --with-debugging=0 --download-suitesparse --download-parmetis
>>>>>>> --download-blacs
>>>>>>> --with-blas-lapack-dir=/opt/intel/system_studio_2015.2.050/mkl
>>>>>>> --CXXFLAGS=-Wl,--no-as-needed --download-scalapack --download-mumps
>>>>>>> --download-metis --download-superlu_dist
>>>>>>> --prefix=/home/dknez/software/libmesh_install/opt_real/petsc
>>>>>>> --download-hypre --download-ml
>>>>>>> [0]PETSC ERROR: #1 PetscKernel_A_gets_inverse_A_5() line 48 in
>>>>>>> /home/dknez/software/petsc-3.6.1/src/mat/impls/baij/seq/dgefa5.c
>>>>>>> [0]PETSC ERROR: #2 MatSOR_SeqAIJ_Inode() line 2808 in
>>>>>>> /home/dknez/software/petsc-3.6.1/src/mat/impls/aij/seq/inode.c
>>>>>>> [0]PETSC ERROR: #3 MatSOR() line 3697 in
>>>>>>> /home/dknez/software/petsc-3.6.1/src/mat/interface/matrix.c
>>>>>>> [0]PETSC ERROR: #4 PCApply_SOR() line 37 in
>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/sor/sor.c
>>>>>>> [0]PETSC ERROR: #5 PCApply() line 482 in
>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c
>>>>>>> [0]PETSC ERROR: #6 KSP_PCApply() line 242 in
>>>>>>> /home/dknez/software/petsc-3.6.1/include/petsc/private/kspimpl.h
>>>>>>> [0]PETSC ERROR: #7 KSPInitialResidual() line 63 in
>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itres.c
>>>>>>> [0]PETSC ERROR: #8 KSPSolve_GMRES() line 235 in
>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/gmres/gmres.c
>>>>>>> [0]PETSC ERROR: #9 KSPSolve() line 604 in
>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
>>>>>>> [0]PETSC ERROR: #10 KSPSolve_Chebyshev() line 381 in
>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/cheby/cheby.c
>>>>>>> [0]PETSC ERROR: #11 KSPSolve() line 604 in
>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
>>>>>>> [0]PETSC ERROR: #12 PCMGMCycle_Private() line 19 in
>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
>>>>>>> [0]PETSC ERROR: #13 PCMGMCycle_Private() line 48 in
>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
>>>>>>> [0]PETSC ERROR: #14 PCApply_MG() line 338 in
>>>>>&

[petsc-users] GAMG and zero pivots follow up

[David Knezevic]

I'm looking into using GAMG, so I wanted to start with a simple 3D
elasticity problem. When I first tried this, I got the following "zero
pivot" error:

---

[0]PETSC ERROR: Zero pivot in LU factorization:
http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot
[0]PETSC ERROR: Zero pivot, row 3
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015
[0]PETSC ERROR: /home/dknez/akselos-dev/scrbe/build/bin/fe_solver-opt_real
on a arch-linux2-c-opt named david-Lenovo by dknez Tue Nov 10 21:26:39 2015
[0]PETSC ERROR: Configure options --with-shared-libraries=1
--with-debugging=0 --download-suitesparse --download-parmetis
--download-blacs
--with-blas-lapack-dir=/opt/intel/system_studio_2015.2.050/mkl
--CXXFLAGS=-Wl,--no-as-needed --download-scalapack --download-mumps
--download-metis --download-superlu_dist
--prefix=/home/dknez/software/libmesh_install/opt_real/petsc
--download-hypre --download-ml
[0]PETSC ERROR: #1 PetscKernel_A_gets_inverse_A_5() line 48 in
/home/dknez/software/petsc-3.6.1/src/mat/impls/baij/seq/dgefa5.c
[0]PETSC ERROR: #2 MatSOR_SeqAIJ_Inode() line 2808 in
/home/dknez/software/petsc-3.6.1/src/mat/impls/aij/seq/inode.c
[0]PETSC ERROR: #3 MatSOR() line 3697 in
/home/dknez/software/petsc-3.6.1/src/mat/interface/matrix.c
[0]PETSC ERROR: #4 PCApply_SOR() line 37 in
/home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/sor/sor.c
[0]PETSC ERROR: #5 PCApply() line 482 in
/home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c
[0]PETSC ERROR: #6 KSP_PCApply() line 242 in
/home/dknez/software/petsc-3.6.1/include/petsc/private/kspimpl.h
[0]PETSC ERROR: #7 KSPInitialResidual() line 63 in
/home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itres.c
[0]PETSC ERROR: #8 KSPSolve_GMRES() line 235 in
/home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/gmres/gmres.c
[0]PETSC ERROR: #9 KSPSolve() line 604 in
/home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #10 KSPSolve_Chebyshev() line 381 in
/home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/cheby/cheby.c
[0]PETSC ERROR: #11 KSPSolve() line 604 in
/home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #12 PCMGMCycle_Private() line 19 in
/home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
[0]PETSC ERROR: #13 PCMGMCycle_Private() line 48 in
/home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
[0]PETSC ERROR: #14 PCApply_MG() line 338 in
/home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
[0]PETSC ERROR: #15 PCApply() line 482 in
/home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c
[0]PETSC ERROR: #16 KSP_PCApply() line 242 in
/home/dknez/software/petsc-3.6.1/include/petsc/private/kspimpl.h
[0]PETSC ERROR: #17 KSPSolve_CG() line 139 in
/home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/cg/cg.c
[0]PETSC ERROR: #18 KSPSolve() line 604 in
/home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c

---

I saw that there was a thread about this in September (subject: "gamg and
zero pivots"), and that the fix is to use "-mg_levels_pc_type jacobi." When
I do that, the solve succeeds (I pasted the -ksp_view at the end of this
email).

So I have two questions about this:

1. Is it surprising that I hit this issue for a 3D elasticity problem? Note
that matrix assembly was done in libMesh, I can look into the structure of
the assembled matrix more carefully, if needed. Also, note that I can solve
this problem with direct solvers just fine.

2. Is there a way to set "-mg_levels_pc_type jacobi" programmatically,
rather than via the command line?

Thanks,
David

---

ksp_view output:


KSP Object: 1 MPI processes type: cg maximum iterations=5000 tolerances:
relative=1e-12, absolute=1e-50, divergence=1 left preconditioning using
nonzero initial guess using PRECONDITIONED norm type for convergence test
PC Object: 1 MPI processes type: gamg MG: type is MULTIPLICATIVE, levels=5
cycles=v Cycles per PCApply=1 Using Galerkin computed coarse grid matrices
GAMG specific options Threshold for dropping small values from graph 0 AGG
specific options Symmetric graph false Coarse grid solver -- level
--- KSP Object: (mg_coarse_) 1 MPI processes
type: gmres GMRES: restart=30, using Classical (unmodified) Gram-Schmidt
Orthogonalization with no iterative refinement GMRES: happy breakdown
tolerance 1e-30 maximum iterations=1, initial guess is zero tolerances:
relative=1e-05, absolute=1e-50, divergence=1 left preconditioning using
NONE norm type for convergence test PC Object: (mg_coarse_) 1 MPI processes
type: bjacobi block Jacobi: number of blocks = 1 Local solve is same for
all blocks, in the following KSP and PC objects: KSP Object:
(mg_coarse_sub_) 1 MPI 

Re: [petsc-users] GAMG and zero pivots follow up

[David Knezevic]

On Tue, Nov 10, 2015 at 10:24 PM, Matthew Knepley <knep...@gmail.com> wrote:

> On Tue, Nov 10, 2015 at 9:21 PM, David Knezevic <
> david.kneze...@akselos.com> wrote:
>
>> On Tue, Nov 10, 2015 at 10:00 PM, Matthew Knepley <knep...@gmail.com>
>> wrote:
>>
>>> On Tue, Nov 10, 2015 at 8:39 PM, David Knezevic <
>>> david.kneze...@akselos.com> wrote:
>>>
>>>> I'm looking into using GAMG, so I wanted to start with a simple 3D
>>>> elasticity problem. When I first tried this, I got the following "zero
>>>> pivot" error:
>>>>
>>>> ---
>>>>
>>>> [0]PETSC ERROR: Zero pivot in LU factorization:
>>>> http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot
>>>> [0]PETSC ERROR: Zero pivot, row 3
>>>> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
>>>> for trouble shooting.
>>>> [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015
>>>> [0]PETSC ERROR:
>>>> /home/dknez/akselos-dev/scrbe/build/bin/fe_solver-opt_real on a
>>>> arch-linux2-c-opt named david-Lenovo by dknez Tue Nov 10 21:26:39 2015
>>>> [0]PETSC ERROR: Configure options --with-shared-libraries=1
>>>> --with-debugging=0 --download-suitesparse --download-parmetis
>>>> --download-blacs
>>>> --with-blas-lapack-dir=/opt/intel/system_studio_2015.2.050/mkl
>>>> --CXXFLAGS=-Wl,--no-as-needed --download-scalapack --download-mumps
>>>> --download-metis --download-superlu_dist
>>>> --prefix=/home/dknez/software/libmesh_install/opt_real/petsc
>>>> --download-hypre --download-ml
>>>> [0]PETSC ERROR: #1 PetscKernel_A_gets_inverse_A_5() line 48 in
>>>> /home/dknez/software/petsc-3.6.1/src/mat/impls/baij/seq/dgefa5.c
>>>> [0]PETSC ERROR: #2 MatSOR_SeqAIJ_Inode() line 2808 in
>>>> /home/dknez/software/petsc-3.6.1/src/mat/impls/aij/seq/inode.c
>>>> [0]PETSC ERROR: #3 MatSOR() line 3697 in
>>>> /home/dknez/software/petsc-3.6.1/src/mat/interface/matrix.c
>>>> [0]PETSC ERROR: #4 PCApply_SOR() line 37 in
>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/sor/sor.c
>>>> [0]PETSC ERROR: #5 PCApply() line 482 in
>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c
>>>> [0]PETSC ERROR: #6 KSP_PCApply() line 242 in
>>>> /home/dknez/software/petsc-3.6.1/include/petsc/private/kspimpl.h
>>>> [0]PETSC ERROR: #7 KSPInitialResidual() line 63 in
>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itres.c
>>>> [0]PETSC ERROR: #8 KSPSolve_GMRES() line 235 in
>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/gmres/gmres.c
>>>> [0]PETSC ERROR: #9 KSPSolve() line 604 in
>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
>>>> [0]PETSC ERROR: #10 KSPSolve_Chebyshev() line 381 in
>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/cheby/cheby.c
>>>> [0]PETSC ERROR: #11 KSPSolve() line 604 in
>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
>>>> [0]PETSC ERROR: #12 PCMGMCycle_Private() line 19 in
>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
>>>> [0]PETSC ERROR: #13 PCMGMCycle_Private() line 48 in
>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
>>>> [0]PETSC ERROR: #14 PCApply_MG() line 338 in
>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
>>>> [0]PETSC ERROR: #15 PCApply() line 482 in
>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c
>>>> [0]PETSC ERROR: #16 KSP_PCApply() line 242 in
>>>> /home/dknez/software/petsc-3.6.1/include/petsc/private/kspimpl.h
>>>> [0]PETSC ERROR: #17 KSPSolve_CG() line 139 in
>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/cg/cg.c
>>>> [0]PETSC ERROR: #18 KSPSolve() line 604 in
>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
>>>>
>>>> ---
>>>>
>>>> I saw that there was a thread about this in September (subject: "gamg
>>>> and zero pivots"), and that the fix is to use "-mg_levels_pc_type
>>>> jacobi." When I do that, the solve succeeds (I pasted the -ksp_view at the
>>>> end of this email).
>>>>
>>>> So I have two questions about

Re: [petsc-users] GAMG and zero pivots follow up

[David Knezevic]

On Tue, Nov 10, 2015 at 10:00 PM, Matthew Knepley <knep...@gmail.com> wrote:

> On Tue, Nov 10, 2015 at 8:39 PM, David Knezevic <
> david.kneze...@akselos.com> wrote:
>
>> I'm looking into using GAMG, so I wanted to start with a simple 3D
>> elasticity problem. When I first tried this, I got the following "zero
>> pivot" error:
>>
>> ---
>>
>> [0]PETSC ERROR: Zero pivot in LU factorization:
>> http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot
>> [0]PETSC ERROR: Zero pivot, row 3
>> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
>> for trouble shooting.
>> [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015
>> [0]PETSC ERROR:
>> /home/dknez/akselos-dev/scrbe/build/bin/fe_solver-opt_real on a
>> arch-linux2-c-opt named david-Lenovo by dknez Tue Nov 10 21:26:39 2015
>> [0]PETSC ERROR: Configure options --with-shared-libraries=1
>> --with-debugging=0 --download-suitesparse --download-parmetis
>> --download-blacs
>> --with-blas-lapack-dir=/opt/intel/system_studio_2015.2.050/mkl
>> --CXXFLAGS=-Wl,--no-as-needed --download-scalapack --download-mumps
>> --download-metis --download-superlu_dist
>> --prefix=/home/dknez/software/libmesh_install/opt_real/petsc
>> --download-hypre --download-ml
>> [0]PETSC ERROR: #1 PetscKernel_A_gets_inverse_A_5() line 48 in
>> /home/dknez/software/petsc-3.6.1/src/mat/impls/baij/seq/dgefa5.c
>> [0]PETSC ERROR: #2 MatSOR_SeqAIJ_Inode() line 2808 in
>> /home/dknez/software/petsc-3.6.1/src/mat/impls/aij/seq/inode.c
>> [0]PETSC ERROR: #3 MatSOR() line 3697 in
>> /home/dknez/software/petsc-3.6.1/src/mat/interface/matrix.c
>> [0]PETSC ERROR: #4 PCApply_SOR() line 37 in
>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/sor/sor.c
>> [0]PETSC ERROR: #5 PCApply() line 482 in
>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c
>> [0]PETSC ERROR: #6 KSP_PCApply() line 242 in
>> /home/dknez/software/petsc-3.6.1/include/petsc/private/kspimpl.h
>> [0]PETSC ERROR: #7 KSPInitialResidual() line 63 in
>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itres.c
>> [0]PETSC ERROR: #8 KSPSolve_GMRES() line 235 in
>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/gmres/gmres.c
>> [0]PETSC ERROR: #9 KSPSolve() line 604 in
>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
>> [0]PETSC ERROR: #10 KSPSolve_Chebyshev() line 381 in
>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/cheby/cheby.c
>> [0]PETSC ERROR: #11 KSPSolve() line 604 in
>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
>> [0]PETSC ERROR: #12 PCMGMCycle_Private() line 19 in
>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
>> [0]PETSC ERROR: #13 PCMGMCycle_Private() line 48 in
>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
>> [0]PETSC ERROR: #14 PCApply_MG() line 338 in
>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c
>> [0]PETSC ERROR: #15 PCApply() line 482 in
>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c
>> [0]PETSC ERROR: #16 KSP_PCApply() line 242 in
>> /home/dknez/software/petsc-3.6.1/include/petsc/private/kspimpl.h
>> [0]PETSC ERROR: #17 KSPSolve_CG() line 139 in
>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/cg/cg.c
>> [0]PETSC ERROR: #18 KSPSolve() line 604 in
>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
>>
>> ---
>>
>> I saw that there was a thread about this in September (subject: "gamg and
>> zero pivots"), and that the fix is to use "-mg_levels_pc_type jacobi."
>> When I do that, the solve succeeds (I pasted the -ksp_view at the end of
>> this email).
>>
>> So I have two questions about this:
>>
>> 1. Is it surprising that I hit this issue for a 3D elasticity problem?
>> Note that matrix assembly was done in libMesh, I can look into the
>> structure of the assembled matrix more carefully, if needed. Also, note
>> that I can solve this problem with direct solvers just fine.
>>
>
> Yes, this seems like a bug, but it could be some strange BC thing I do not
> understand.
>


OK, I can look into the matrix in more detail. I agree that it should have
a non-zero diagonal, so I'll have a look at what's happening with that.




> Naively, the elastic element matrix has a nonzero diagonal. I see that you
> are doing LU
> of size 5. That seems strange for 3D elasticity. Am I missing something? I
> would expec

[petsc-users] Automatically re-solving after MUMPS error

[David Knezevic]

In some cases, I get MUMPS error -9, i.e.:
[2]PETSC ERROR: Error reported by MUMPS in numerical factorization phase:
INFO(1)=-9, INFO(2)=98927

This is easily fixed by re-running the executable with -mat_mumps_icntl_14
on the commandline.

However, I would like to update my code in order to do this automatically,
i.e. detect the -9 error and re-run with the appropriate option. Is there a
recommended way to do this? It seems to me that I could do this with a
PETSc error handler (e.g. PetscPushErrorHandler) in order to call a
function that sets the appropriate option and solves again, is that right?
Are there any examples that illustrate this type of thing?

Thanks,
David

 

Re: [petsc-users] Automatically re-solving after MUMPS error

[David Knezevic]

On Tue, Sep 15, 2015 at 7:29 PM, Matthew Knepley <knep...@gmail.com> wrote:

> On Tue, Sep 15, 2015 at 4:30 AM, David Knezevic <
> david.kneze...@akselos.com> wrote:
>
>> In some cases, I get MUMPS error -9, i.e.:
>> [2]PETSC ERROR: Error reported by MUMPS in numerical factorization phase:
>> INFO(1)=-9, INFO(2)=98927
>>
>> This is easily fixed by re-running the executable with
>> -mat_mumps_icntl_14 on the commandline.
>>
>> However, I would like to update my code in order to do this
>> automatically, i.e. detect the -9 error and re-run with the appropriate
>> option. Is there a recommended way to do this? It seems to me that I could
>> do this with a PETSc error handler (e.g. PetscPushErrorHandler) in order to
>> call a function that sets the appropriate option and solves again, is that
>> right? Are there any examples that illustrate this type of thing?
>>
>
> I would not use the error handler. I would just check the ierr return code
> from the solver. I think you need the
> INFO output, for which you can use MatMumpsGetInfo().
>


OK, that sounds good (and much simpler than what I had in mind), thanks for
the help!

David

[petsc-users] Adding rows and columns to a matrix

[David Knezevic]

I need to add Lagrange multipliers to an existing system of equations,
which requires me to expand the matrix.

I figured I would allocate a new larger matrix with the appropriate
non-zero structure, and then copy the entries from the smaller matrix over
to it. But I was wondering if there might be a better way to do this?

Thanks,
David

Re: [petsc-users] Adding rows and columns to a matrix

[David Knezevic]

On Thu, Sep 3, 2015 at 8:11 PM, Barry Smith <bsm...@mcs.anl.gov> wrote:

>
>   You might also consider the MatNest object. It's job is to help
> construct matrices that have a natural set of sub matrices such as yours.
> Also look at MatGetLocalSubMat()
>


OK, thanks for the pointers.

Best regards,
David




>
> > On Sep 3, 2015, at 6:52 PM, David Knezevic <david.kneze...@akselos.com>
> wrote:
> >
> > I need to add Lagrange multipliers to an existing system of equations,
> which requires me to expand the matrix.
> >
> > I figured I would allocate a new larger matrix with the appropriate
> non-zero structure, and then copy the entries from the smaller matrix over
> to it. But I was wondering if there might be a better way to do this?
> >
> > Thanks,
> > David
> >
> >
>
>

Re: [petsc-users] Variatonal inequalities

[David Knezevic]

On Sat, Aug 22, 2015 at 10:29 PM, Barry Smith bsm...@mcs.anl.gov wrote:


   David,

Currently the only way to do this without adding a lot of additional
 PETSc code is to add additional variables such that only box constraints
 appear in the final problem. For example say you have constraints   c = Ax
 = d then introduce new variables y = Ax and then you have the larger
 problem of unknowns (x,y) and box constrains on y with -infinity and
 +infinity constraints on x.



OK, that makes sense, thanks for the info!

David

[petsc-users] Variatonal inequalities

[David Knezevic]

Hi all,

I see from Section 5.7 of the manual that SNES supports box constraints on
variables, which is great. However, I was also hoping to also be able to
consider general linear inequality constraints, so I was wondering if
anyone has any suggestions on how (or if) that could be done with PETSc?

Thanks,
David

Re: [petsc-users] MUMPS error

[David Knezevic]

On Sat, May 16, 2015 at 8:08 AM, venkatesh g venkateshg...@gmail.com
wrote:

 Hi,
 I am trying to solving AX=lambda BX eigenvalue problem.

 A and B are of sizes 3600x3600

 I run with this command :

 'mpiexec -np 4 ./ex7 -f1 a2 -f2 b2 -eps_nev 1 -st_type sinvert -eps_max_it
 5000 -st_ksp_type preonly -st_pc_type lu -st_pc_factor_mat_solver_package
 mumps'

 I get this error: (I get result only when I give 1 or 2 processors)
 Reading COMPLEX matrices from binary files...
 [0]PETSC ERROR: - Error Message
 
 [0]PETSC ERROR: Error in external library!
 [0]PETSC ERROR: Error reported by MUMPS in numerical factorization phase:
 INFO(1)=-9, INFO(2)=2024



The MUMPS error types are described in Chapter 7 of the MUMPS manual. In
this case you have INFO(1)=-9, which is explained in the manual as:

–9 Main internal real/complex workarray S too small. If INFO(2) is
positive, then the number of entries that are missing in S at the moment
when the error is raised is available in INFO(2). If INFO(2) is negative,
then its absolute value should be multiplied by 1 million. If an error –9
occurs, the user should increase the value of ICNTL(14) before calling the
factorization (JOB=2) again, except if ICNTL(23) is provided, in which case
ICNTL(23) should be increased.

This says that you should use ICTNL(14) to increase the working space size:

ICNTL(14) is accessed by the host both during the analysis and the
factorization phases. It corresponds to the percentage increase in the
estimated working space. When significant extra fill-in is caused by
numerical pivoting, increasing ICNTL(14) may help. Except in special cases,
the default value is 20 (which corresponds to a 20 % increase).

So, for example, you can avoid this error via the following command line
argument to PETSc: -mat_mumps_icntl_14 30, where 30 indicates that we
allow a 30% increase in the workspace instead of the default 20%.

David

Re: [petsc-users] TAO IPM

[David Knezevic]

On Thu, Apr 23, 2015 at 2:57 PM, Jason Sarich jason.sar...@gmail.com
wrote:

 Hi David,

 The IPM code is a naive implementation and not very battle-tested. There
 is nothing specifically broken, but the linear solve operator becomes very
 ill-conditioned and for most cases will require a direct solver (e.g.
 superlu). This is the major problem, and if you're ok with that then it may
 work for you.
 Unfortunately, we don't have anything else to recommend for generally
 constrained problems.



A quick follow-up question: I notice that the toy example I referred to
earlier in this thread uses VecCreateSeq and MatCreateSeq, so I just wanted
to check if TAO IPM is restricted to sequential data? (My guess is that
it's not restricted to sequential data, and that Seq is used only because
the problem is so small, but I just wanted to check.)

Thanks,
David

Re: [petsc-users] TAO IPM

[David Knezevic]

On Thu, Apr 23, 2015 at 2:57 PM, Jason Sarich jason.sar...@gmail.com
wrote:

 Hi David,

 The IPM code is a naive implementation and not very battle-tested. There
 is nothing specifically broken, but the linear solve operator becomes very
 ill-conditioned and for most cases will require a direct solver (e.g.
 superlu). This is the major problem, and if you're ok with that then it may
 work for you.
 Unfortunately, we don't have anything else to recommend for generally
 constrained problems.



Hi Jason,

OK, understood, thanks very much for that info.

Best regards,
David

[petsc-users] TAO IPM

[David Knezevic]

I'm interested in trying out TAO IPM for constrained optimization with
inequality, equality and bound constraints. I implemented a test case based
on the toy example
http://www.mcs.anl.gov/petsc/petsc-current/src/tao/constrained/examples/tutorials/toy.c.html
and it works well.

But I noticed the following comment in the documentation
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Tao/TAOIPM.html
:
This algorithm is more of a place-holder for future constrained
optimization algorithms and should not yet be used for large problems or
production code.

So I was wondering what the status of this is? Is there anything that
specifically doesn't work, or is it just not production-ready at this
stage? (I didn't see any documentation about IPM in the TAO manual, so I
figured I'd ask here).

If IPM isn't recommended at this stage, is there anything else in TAO that
works for generally constrained problems?

Thanks,
David

Re: [petsc-users] Exploiting symmetry with direct solvers

[David Knezevic]

On Wed, Apr 8, 2015 at 9:28 PM, Hong hzh...@mcs.anl.gov wrote:

 David ,

 I'm curious about how to do a symmetric LDL^T factorization (instead of
 LU) with MUMPS and SuperLU. Based on this example:


 http://www.mcs.anl.gov/petsc/petsc-3.4/src/ksp/ksp/examples/tutorials/ex52.c.html

 my understanding is as follows:

 - With MUMPS I gather that we need to specify:

 MatSetOption(A,MAT_SPD,PETSC_TRUE);
 PCSetType(pc,PCCHOLESKY);

 I guess -pc_type cholesky on the command line is equivalent to the
 PCSetType call, right? Is specifying MAT_SPD required in order for MUMPS to
 do an LDL^T factorization?

 Mumps supports Cholesky factorization for symmetric, and symmetric+spd
 matrices.  You may consult mumps user manual.



OK, thanks. But I was also wondering about the PETSc interface to MUMPS.

I consulted the MUMPS manual, but it just says that you need to specify
SYM=1 in order to get an LDL^T factorization. I gather that if we set
MatSetOption(A,MAT_SPD,PETSC_TRUE); in PETSc then that will tell MUMPS to
use SYM=1, right?

Thanks,
David

Re: [petsc-users] Exploiting symmetry with direct solvers

[David Knezevic]

On Wed, Apr 8, 2015 at 11:06 PM, Hong hzh...@mcs.anl.gov wrote:

 David:


 Mumps supports Cholesky factorization for symmetric, and symmetric+spd
 matrices.  You may consult mumps user manual.



 OK, thanks. But I was also wondering about the PETSc interface to MUMPS.

 I consulted the MUMPS manual, but it just says that you need to specify
 SYM=1 in order to get an LDL^T factorization. I gather that if we set
 MatSetOption(A,MAT_SPD,PETSC_TRUE); in PETSc then that will tell MUMPS to
 use SYM=1, right?


 When user requests Cholesky factorization, PETSc sets  SYM=2 as default;
 If user sets flag of spd, then the interface uses SYM=1 in order to get
 an LDL^T factorization. See MatGetFactor_xxx_mumps() in
 petsc/src/mat/impls/aij/mpi/mumps/mumps.c:

 ...
  B-factortype = MAT_FACTOR_CHOLESKY;
   if (A-spd_set  A-spd) mumps-sym = 1;
   else  mumps-sym = 2;
 ...



OK, got it, thanks!

David

[petsc-users] Exploiting symmetry with direct solvers

[David Knezevic]

I'm curious about how to do a symmetric LDL^T factorization (instead of LU)
with MUMPS and SuperLU. Based on this example:

http://www.mcs.anl.gov/petsc/petsc-3.4/src/ksp/ksp/examples/tutorials/ex52.c.html

my understanding is as follows:

- With MUMPS I gather that we need to specify:

MatSetOption(A,MAT_SPD,PETSC_TRUE);
PCSetType(pc,PCCHOLESKY);

I guess -pc_type cholesky on the command line is equivalent to the
PCSetType call, right? Is specifying MAT_SPD required in order for MUMPS to
do an LDL^T factorization?

- I gather that SuperLU doesn't provide a symmetric factorization.

Is this correct?

Thanks!
David

[petsc-users] Status of mkl_cpardiso configuration option?

[David Knezevic]

Hi all,

I'm interested in trying out mkl_cpardiso with PETSc (since I gather that
mkl_pardiso doesn't work with MPI, and mkl_cpardiso is included in MKL). I
saw a thread on the petsc-users list about this from 2013 (subject:
mkl-cpardiso configuration), so I was wondering what the status of this
is?

I see that there is a configuration option for mkl_pardiso (I tried this
and it works well for me). But I gather that there is no configuration
option for mkl_cpardiso at this stage, is that correct? Is this
configuration option something that you're interested in adding? If so I'd
be happy to help out if possible.

Thanks,
David

Re: [petsc-users] Status of mkl_cpardiso configuration option?

[David Knezevic]

On Mon, Apr 6, 2015 at 12:19 PM, Satish Balay ba...@mcs.anl.gov wrote:

 If you use petsc development branch 'master' from git - you should be
 able to use mkl cPardiso.


 https://bitbucket.org/petsc/petsc/commits/d305a81be0fbb2fe8c02e3da7510a1143eafa263

 Satish



OK, thanks!

David

Re: [petsc-users] Error configuring with --download-ml and Intel MKL

[David Knezevic]

On Thu, Apr 2, 2015 at 1:52 AM, Satish Balay ba...@mcs.anl.gov wrote:

 And to get g++ to behave similar to gcc on 14.04 - you can add an
 extra configure option LIBS=-Wl,--no-as-needed -lc [so
 --with-clanguage=cxx should work with it]

 ./configure --with-scalar-type=complex --with-clanguage=cxx
 --download-metis --download-parmetis --download-superlu_dist
 --with-blas-lapack-dir=/opt/intel/system_studio_2015.2.050/mkl/lib/intel64
 LIBS=-Wl,--no-as-needed -lc



Thanks! Adding the LIBS option fixed this for me.

David

Re: [petsc-users] Error configuring with --download-ml and Intel MKL

[David Knezevic]

On Thu, Apr 2, 2015 at 11:58 AM, Satish Balay ba...@mcs.anl.gov wrote:

 On Thu, 2 Apr 2015, David Knezevic wrote:

  On Thu, Apr 2, 2015 at 1:52 AM, Satish Balay ba...@mcs.anl.gov wrote:
 
   And to get g++ to behave similar to gcc on 14.04 - you can add an
   extra configure option LIBS=-Wl,--no-as-needed -lc [so
   --with-clanguage=cxx should work with it]
  
   ./configure --with-scalar-type=complex --with-clanguage=cxx
   --download-metis --download-parmetis --download-superlu_dist
  
 --with-blas-lapack-dir=/opt/intel/system_studio_2015.2.050/mkl/lib/intel64
   LIBS=-Wl,--no-as-needed -lc
  
 
 
  Thanks! Adding the LIBS option fixed this for me.

 Looks like you can also use:

 CXXFLAGS=-Wl,--no-as-needed

 And --download-ml should work with this.. [as CXXFLAGS get passed down to
 ML's configure]



Yep, that works for me for both ML and SuperLU with MKL. Thanks so much for
your help!

David

Re: [petsc-users] Error configuring with --download-ml and Intel MKL

[David Knezevic]

On Thu, Mar 26, 2015 at 6:23 PM, Matthew Knepley knep...@gmail.com wrote:

 On Thu, Mar 26, 2015 at 5:10 PM, David Knezevic 
 david.kneze...@akselos.com wrote:

 On Thu, Mar 26, 2015 at 4:21 PM, Matthew Knepley knep...@gmail.com
 wrote:

 On Thu, Mar 26, 2015 at 1:48 PM, David Knezevic 
 david.kneze...@akselos.com wrote:

 Hi all,

 I'm trying to configure PETSc using Intel's MKL and with --download-ml.
 Here is my configure line:

 ./configure
 --with-blas-lapack-dir=/opt/intel/composer_xe_2015/mkl/lib/intel64
 --download-ml

 I get an error when ML does checking for dgemm. The configure.log is
 attached. I was wondering if anyone has any suggestions about how I can get
 this to work?


 We will need $PETSC_ARCH/externalpackages/ml-*/config.log to see where
 their test failed.



 I've attached this file.


 I do not understand the error:

 configure:10865: mpicxx -o conftest  -Wall -Wwrite-strings
 -Wno-strict-aliasing -Wno-unknown-pragmas -g -O0   -fPIC
 -DMPICH_SKIP_MPICXX -DOMPI_SKIP_MPICXX -I/usr/lib/openmpi/include
 -I/usr/lib/openmpi/include/openmpi -g -O2 -I/usr/lib/openmpi/include
 -I/usr/lib/openmpi/include/openmpi  conftest.cpp
 -Wl,-rpath,/opt/intel/composer_xe_2015/mkl/lib/intel64
 -L/opt/intel/composer_xe_2015/mkl/lib/intel64 -lmkl_intel_lp64
 -lmkl_sequential -lmkl_core -lpthread -lm -Wl,-rpath,/usr/lib/openmpi/lib
 -L/usr/lib/openmpi/lib -Wl,-rpath,/usr/lib/gcc/x86_64-linux-gnu/4.8
 -L/usr/lib/gcc/x86_64-linux-gnu/4.8 -Wl,-rpath,/usr/lib/x86_64-linux-gnu
 -L/usr/lib/x86_64-linux-gnu -Wl,-rpath,/lib/x86_64-linux-gnu
 -L/lib/x86_64-linux-gnu -lmpi_f90 -lmpi_f77 -lgfortran -lm
 -Wl,-rpath,/usr/lib/openmpi/lib
 -Wl,-rpath,/usr/lib/gcc/x86_64-linux-gnu/4.8
 -Wl,-rpath,/usr/lib/x86_64-linux-gnu -Wl,-rpath,/lib/x86_64-linux-gnu
 -lgfortran -lm -lquadmath -lm   -L/usr//lib -L/usr/lib/openmpi/lib
 -L/usr/lib/gcc/x86_64-linux-gnu/4.8
 -L/usr/lib/gcc/x86_64-linux-gnu/4.8/../../../x86_64-linux-gnu
 -L/usr/lib/gcc/x86_64-linux-gnu/4.8/../../../../lib -L/lib/x86_64-linux-gnu
 -L/lib/../lib -L/usr/lib/x86_64-linux-gnu -L/usr/lib/../lib
 -L/usr/lib/gcc/x86_64-linux-gnu/4.8/../../.. -lmpi_f90 -lmpi_f77 -lmpi -ldl
 -lhwloc -lgfortran -lm -lquadmath -lpthread 5
 /opt/intel/composer_xe_2015/mkl/lib/intel64/libmkl_core.so: undefined
 reference to `logf'
 /opt/intel/composer_xe_2015/mkl/lib/intel64/libmkl_core.so: undefined
 reference to `atan2'
 /opt/intel/composer_xe_2015/mkl/lib/intel64/libmkl_core.so: undefined
 reference to `sin'
 /opt/intel/composer_xe_2015/mkl/lib/intel64/libmkl_core.so: undefined
 reference to `fabs'
 /opt/intel/composer_xe_2015/mkl/lib/intel64/libmkl_core.so: undefined
 reference to `exp'
 /opt/intel/composer_xe_2015/mkl/lib/intel64/libmkl_core.so: undefined
 reference to `cos'
 /opt/intel/composer_xe_2015/mkl/lib/intel64/libmkl_core.so: undefined
 reference to `sqrt'
 /opt/intel/composer_xe_2015/mkl/lib/intel64/libmkl_intel_lp64.so:
 undefined reference to `log'
 /opt/intel/composer_xe_2015/mkl/lib/intel64/libmkl_core.so: undefined
 reference to `pow'
 /opt/intel/composer_xe_2015/mkl/lib/intel64/libmkl_core.so: undefined
 reference to `log10'
 /opt/intel/composer_xe_2015/mkl/lib/intel64/libmkl_core.so: undefined
 reference to `ceil'
 /opt/intel/composer_xe_2015/mkl/lib/intel64/libmkl_core.so: undefined
 reference to `expf'
 collect2: error: ld returned 1 exit status

 Clearly these symbols should be in -lm, which is in the link line, and it
 linked when we tried in our configure. Can
 you run this link line manually and figure out what is wrong?



I didn't get to the bottom of this issue with ML+MKL yet, but I have
another case that's more pressing right now, and I was wondering if you
might be able to help me out with it. The configuration I'm trying is as
follows:

./configure --with-scalar-type=complex --with-clanguage=cxx
--download-metis --download-parmetis --download-superlu_dist
--with-blas-lapack-dir=/path/to/mkl

So I'm still trying to use MKL for the time being (though I certainly
appreciate your input about the drawbacks of this!), along with
SuperLU_dist in complex mode. This configuration fails for me. I don't see
any configuration log file in $PETSC_ARCH/externalpackages/SuperLU_DIST, is
it supposed to be there somewhere?

The error in the PETSc configure.log is pasted below. Interestingly if I
configure in real mode (by removing --with-scalar-type=complex
--with-clanguage=cxx) then it works fine.

If anyone has any suggestions about what I might be able to do to fix this,
that'd be most appreciated.

Thanks!
David


***
 UNABLE to CONFIGURE with GIVEN OPTIONS(see configure.log for
details):
---
Downloaded superlu_dist could not be used. Please check install in
/home/dknez/software/petsc-3.5.2/arch-linux2-cxx-debug
***
  File

Re: [petsc-users] Error configuring with --download-ml and Intel MKL

[David Knezevic]

On Thu, Mar 26, 2015 at 4:21 PM, Matthew Knepley knep...@gmail.com wrote:

 On Thu, Mar 26, 2015 at 1:48 PM, David Knezevic 
 david.kneze...@akselos.com wrote:

 Hi all,

 I'm trying to configure PETSc using Intel's MKL and with --download-ml.
 Here is my configure line:

 ./configure
 --with-blas-lapack-dir=/opt/intel/composer_xe_2015/mkl/lib/intel64
 --download-ml

 I get an error when ML does checking for dgemm. The configure.log is
 attached. I was wondering if anyone has any suggestions about how I can get
 this to work?


 We will need $PETSC_ARCH/externalpackages/ml-*/config.log to see where
 their test failed.



I've attached this file.




 However, MKL is almost never worth it.


OK, interesting, why not? The reason I'm using MKL in this case is because
I found that the direct solvers MUMPS and SuperLU were surprising slow on a
Haswell server. After rebuilding PETSc with MKL, these solvers were more
than 3x faster in some cases.

David
This file contains any messages produced by compilers while
running configure, to aid debugging if configure makes a mistake.

It was created by ml configure 6.2, which was
generated by GNU Autoconf 2.61.  Invocation command line was

  $ ./configure --prefix=/home/akselos/software/petsc-3.5.3/arch-linux2-c-debug --libdir=/home/akselos/software/petsc-3.5.3/arch-linux2-c-debug/lib --disable-ml-epetra --disable-ml-aztecoo --disable-ml-examples --disable-tests CC=mpicc --with-cflags= -fPIC -Wall -Wwrite-strings -Wno-strict-aliasing -Wno-unknown-pragmas -g3 -O0  -DMPICH_SKIP_MPICXX -DOMPI_SKIP_MPICXX -I/usr/lib/openmpi/include -I/usr/lib/openmpi/include/openmpi CPPFLAGS=-I/usr/lib/openmpi/include -I/usr/lib/openmpi/include/openmpi F77=mpif90 --with-fflags= -fPIC -Wall -Wno-unused-variable -ffree-line-length-0 -Wno-unused-dummy-argument -g -O0   -I/usr/lib/openmpi/include -I/usr/lib/openmpi/include/openmpi CXX=mpicxx --with-cxxflags= -Wall -Wwrite-strings -Wno-strict-aliasing -Wno-unknown-pragmas -g -O0   -fPIC   -DMPICH_SKIP_MPICXX -DOMPI_SKIP_MPICXX -I/usr/lib/openmpi/include -I/usr/lib/openmpi/include/openmpi --enable-mpi --with-mpi-libs= --with-blas=-Wl,-rpath,/opt/intel/composer_xe_2015/mkl/lib/intel64 -L/opt/intel/composer_xe_2015/mkl/lib/intel64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -Wl,-rpath,/usr/lib/openmpi/lib -L/usr/lib/openmpi/lib -Wl,-rpath,/usr/lib/gcc/x86_64-linux-gnu/4.8 -L/usr/lib/gcc/x86_64-linux-gnu/4.8 -Wl,-rpath,/usr/lib/x86_64-linux-gnu -L/usr/lib/x86_64-linux-gnu -Wl,-rpath,/lib/x86_64-linux-gnu -L/lib/x86_64-linux-gnu -lmpi_f90 -lmpi_f77 -lgfortran -lm -Wl,-rpath,/usr/lib/openmpi/lib -Wl,-rpath,/usr/lib/gcc/x86_64-linux-gnu/4.8 -Wl,-rpath,/usr/lib/x86_64-linux-gnu -Wl,-rpath,/lib/x86_64-linux-gnu -lgfortran -lm -lquadmath -lm

## - ##
## Platform. ##
## - ##

hostname = build-scrbe
uname -m = x86_64
uname -r = 3.16.0-31-generic
uname -s = Linux
uname -v = #43~14.04.1-Ubuntu SMP Tue Mar 10 20:13:38 UTC 2015

/usr/bin/uname -p = unknown
/bin/uname -X = unknown

/bin/arch  = unknown
/usr/bin/arch -k   = unknown
/usr/convex/getsysinfo = unknown
/usr/bin/hostinfo  = unknown
/bin/machine   = unknown
/usr/bin/oslevel   = unknown
/bin/universe  = unknown

PATH: /usr/local/sbin
PATH: /usr/local/bin
PATH: /usr/sbin
PATH: /usr/bin
PATH: /sbin
PATH: /bin
PATH: /usr/games
PATH: /usr/local/games


## --- ##
## Core tests. ##
## --- ##

configure:2160: checking whether to enable maintainer-specific portions of Makefiles
configure:2169: result: no
configure:2191: checking build system type
configure:2209: result: x86_64-unknown-linux-gnu
configure:2231: checking host system type
configure:2246: result: x86_64-unknown-linux-gnu
configure:2268: checking target system type
configure:2283: result: x86_64-unknown-linux-gnu
configure:2329: checking for a BSD-compatible install
configure:2385: result: /usr/bin/install -c
configure:2396: checking whether build environment is sane
configure:2439: result: yes
configure:2467: checking for a thread-safe mkdir -p
configure:2506: result: /bin/mkdir -p
configure:2519: checking for gawk
configure:2535: found /usr/bin/gawk
configure:2546: result: gawk
configure:2557: checking whether make sets $(MAKE)
configure:2578: result: yes
configure:2752: checking how to create a ustar tar archive
configure:2765: tar --version
tar (GNU tar) 1.27.1
Copyright (C) 2013 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later http://gnu.org/licenses/gpl.html.
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.

Written by John Gilmore and Jay Fenlason.
configure:2768: $? = 0
configure:2808: tardir=conftest.dir  eval tar --format=ustar -chf - $tardir conftest.tar
configure:2811: $? = 0
configure:2815: tar -xf - conftest.tar
configure:2818: $? = 0
configure:2831: result: gnutar
configure:2993: checking user-defined MPI libraries
configure:2995: result: 
configure

Re: [petsc-users] Near nullspace for preconditioners other than GAMG?

[David Knezevic]

On Tue, Mar 10, 2015 at 5:50 PM, Mark Adams mfad...@lbl.gov wrote:

 Penalty methods are hell on solvers.  There is a good change that GAMG
 will do the right thing with a finite threshold. and keep this stiff
 section separate.  You might want to start wit a small number of levels and
 a large coarse grid and get it working.  Then add more levels and see if it
 breaks.



OK, thanks for the help, much appreciated. I'll look into this.

Best regards,
David








 On Sun, Mar 8, 2015 at 11:38 AM, David Knezevic 
 david.kneze...@akselos.com wrote:

 On Sun, Mar 8, 2015 at 10:31 AM, Mark Adams mfad...@lbl.gov wrote:




 Another situation I've been meaning to ask about: If the elasticity
 model includes some near rigid regions (very high Young's modulus)


 Do you mean Poisson ratio?



 No, just high Young's modulus in subregions. Engineers sometimes
 implement a rotation boundary condition about a node by including rigid
 (or, in practice, very high Young's modulus) parts in a model. For example,
 imagine including a rigid pyramid in a model, in which the tip of the
 pyramid is is clamped so that the base of the pyramid (which is connected
 to a surface in the model) can rotate about the tip.

 This works fine with a direct solver, so I was curious about getting a
 good iterative solver for this case too.



 Nearly incompressible is hard.  bigger smoothers, like ASM, can help.


 Agreed. But the Poisson ratio is still just 0.3 or so in this context, so
 locking and near incompressibility aren't an issue here.

 Thanks,
 David





 On Sat, Mar 7, 2015 at 5:57 PM, Mark Adams mfad...@lbl.gov wrote:

 FYI, stüben used classical AMG for elasticity but he has articulated
 his code for elasticity more than Hypre as I understand it.  Hypre can work
 OK for elasticity in my experience.  Its worth a try.

 Mark

 On Thu, Mar 5, 2015 at 5:27 PM, David Knezevic 
 david.kneze...@akselos.com wrote:

 OK, got it, thanks!

 David


 On Thu, Mar 5, 2015 at 5:08 PM, Jed Brown j...@jedbrown.org wrote:

 David Knezevic david.kneze...@akselos.com writes:
  I was just wondering if its possible to achieve the same sort of
 thing with
  other AMG solvers (e.g. BoomerAMG)? I assume that
 MatSetNearNullSpace does
  nothing for external solvers like hypre, right?

 It is used by ML (smoothed aggregation), but not BoomerAMG (classical
 AMG) which uses an algorithm that doesn't have a natural place for
 such
 information.  To my knowledge, classical AMG is not widely used for
 elasticity.  It is very robust for M-matrices.

Re: [petsc-users] Near nullspace for preconditioners other than GAMG?

[David Knezevic]

On Sun, Mar 8, 2015 at 10:31 AM, Mark Adams mfad...@lbl.gov wrote:




 Another situation I've been meaning to ask about: If the elasticity model
 includes some near rigid regions (very high Young's modulus)


 Do you mean Poisson ratio?



No, just high Young's modulus in subregions. Engineers sometimes implement
a rotation boundary condition about a node by including rigid (or, in
practice, very high Young's modulus) parts in a model. For example, imagine
including a rigid pyramid in a model, in which the tip of the pyramid is is
clamped so that the base of the pyramid (which is connected to a surface in
the model) can rotate about the tip.

This works fine with a direct solver, so I was curious about getting a good
iterative solver for this case too.



 Nearly incompressible is hard.  bigger smoothers, like ASM, can help.


Agreed. But the Poisson ratio is still just 0.3 or so in this context, so
locking and near incompressibility aren't an issue here.

Thanks,
David





 On Sat, Mar 7, 2015 at 5:57 PM, Mark Adams mfad...@lbl.gov wrote:

 FYI, stüben used classical AMG for elasticity but he has articulated his
 code for elasticity more than Hypre as I understand it.  Hypre can work OK
 for elasticity in my experience.  Its worth a try.

 Mark

 On Thu, Mar 5, 2015 at 5:27 PM, David Knezevic 
 david.kneze...@akselos.com wrote:

 OK, got it, thanks!

 David


 On Thu, Mar 5, 2015 at 5:08 PM, Jed Brown j...@jedbrown.org wrote:

 David Knezevic david.kneze...@akselos.com writes:
  I was just wondering if its possible to achieve the same sort of
 thing with
  other AMG solvers (e.g. BoomerAMG)? I assume that MatSetNearNullSpace
 does
  nothing for external solvers like hypre, right?

 It is used by ML (smoothed aggregation), but not BoomerAMG (classical
 AMG) which uses an algorithm that doesn't have a natural place for such
 information.  To my knowledge, classical AMG is not widely used for
 elasticity.  It is very robust for M-matrices.

Re: [petsc-users] Near nullspace for preconditioners other than GAMG?

[David Knezevic]

OK, thanks for letting me know.

I've tried GAMG and ML with MatNullSpaceCreateRigidBody and those both work
well for me. (I tried ML after Jed pointed out that ML also uses the near
nullspace.)

I have also tried hypre but the convergence hasn't been as good for the
elasticity models I've been using, though I may not have been setting the
hypre options in an optimal way for elasticity.

Another situation I've been meaning to ask about: If the elasticity model
includes some near rigid regions (very high Young's modulus) all
iterative solvers I've tried fare poorly. I guess it's because the highly
contrasting stiffnesses give a large condition number. Is there a standard
way to get the preconditioner compensate for this?

Thanks,
David



On Sat, Mar 7, 2015 at 5:57 PM, Mark Adams mfad...@lbl.gov wrote:

 FYI, stüben used classical AMG for elasticity but he has articulated his
 code for elasticity more than Hypre as I understand it.  Hypre can work OK
 for elasticity in my experience.  Its worth a try.

 Mark

 On Thu, Mar 5, 2015 at 5:27 PM, David Knezevic david.kneze...@akselos.com
  wrote:

 OK, got it, thanks!

 David


 On Thu, Mar 5, 2015 at 5:08 PM, Jed Brown j...@jedbrown.org wrote:

 David Knezevic david.kneze...@akselos.com writes:
  I was just wondering if its possible to achieve the same sort of thing
 with
  other AMG solvers (e.g. BoomerAMG)? I assume that MatSetNearNullSpace
 does
  nothing for external solvers like hypre, right?

 It is used by ML (smoothed aggregation), but not BoomerAMG (classical
 AMG) which uses an algorithm that doesn't have a natural place for such
 information.  To my knowledge, classical AMG is not widely used for
 elasticity.  It is very robust for M-matrices.

[petsc-users] Near nullspace for preconditioners other than GAMG?

[David Knezevic]

Hi all,

I found that using MatNullSpaceCreateRigidBody as per example 49
http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex49.c.html
works well with GAMG for a standard elasticity problem, so that's great.

I was just wondering if its possible to achieve the same sort of thing with
other AMG solvers (e.g. BoomerAMG)? I assume that MatSetNearNullSpace does
nothing for external solvers like hypre, right?

Thanks,
David

 http://www.akselos.com

Re: [petsc-users] Near nullspace for preconditioners other than GAMG?

[David Knezevic]

OK, got it, thanks!

David


On Thu, Mar 5, 2015 at 5:08 PM, Jed Brown j...@jedbrown.org wrote:

 David Knezevic david.kneze...@akselos.com writes:
  I was just wondering if its possible to achieve the same sort of thing
 with
  other AMG solvers (e.g. BoomerAMG)? I assume that MatSetNearNullSpace
 does
  nothing for external solvers like hypre, right?

 It is used by ML (smoothed aggregation), but not BoomerAMG (classical
 AMG) which uses an algorithm that doesn't have a natural place for such
 information.  To my knowledge, classical AMG is not widely used for
 elasticity.  It is very robust for M-matrices.

Re: [petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length

[David Knezevic]

Hi Dmitry,

Thanks so much, that's great! I've updated the pull request now to include
your patch:
https://github.com/dknez/libmesh/commit/247cc12f7536a5f848fd5c63765d97af307d9fa7
We can update the code later, once MatReset() is available in PETSc.

Regarding the convergence issues: I think the main issue is that here we're
using the penalty method for contact, which gives a non-smooth term in the
PDE, and hence we don't get nice convergence with the nonlinear solver. In
particular, the SNES solver terminates with DIVERGED_LINE_SEARCH, but I
don't think that's very surprising here. The simulation results seem to be
good though (e.g. see the screenshots in the PR). I'd be interested in your
thoughts on this though?

I'm not sure why the KSP isn't converging fully. I didn't put any thought
into the choice of default solver, happy to change that (to CG or
whatever). But I also used a direct solver for comparison, and SuperLU
gives the same overall result that I get with GMRES (though MUMPS crashes
with Numerically singular matrix for some reason).

Regarding the VI solver: That sounds very good, I'll be interested to try
the overhauled VI solver out for contact problems, once it's available.

Thanks,
David



On Wed, Feb 25, 2015 at 9:56 AM, Dmitry Karpeyev karp...@mcs.anl.gov
wrote:

 The previous patch didn't clean out the Mat data structure thoroughly
 enough. Please try the attached patch.  It runs for me, but there is no
 convergence.  In particular, KSP doesn't seem to converge, although that's
 with the default GMRES+BJACOBI/ILU(0) options. So I don't think resetting
 the Mat is the cause.

 Let me know if this works for you. MatReset() will be in the next release
 as well as, hopefully, an overhauled VI solver that should be able to
 handle these contact problems.

 Dmitry.

 On Mon Feb 23 2015 at 10:17:19 AM David Knezevic 
 david.kneze...@akselos.com wrote:

 OK, sounds good. Let me know if I can help with the digging.

 David



 On Mon, Feb 23, 2015 at 11:10 AM, Dmitry Karpeyev karp...@mcs.anl.gov
 wrote:

 I just tried building against petsc/master, but there needs to be more
 work on libMesh before it can work with petsc/master:
 the new VecLockPush()/Pop() stuff isn't respected by vector manipulation
 in libMesh.
 I put a hack equivalent to MatReset() into your branch (patch attached,
 in case you want to take a look at it),
 but it generates the same error in MatCreateColmap that you reported
 earlier.  It's odd that it occurs
 on the second nonlinear iteration.   I'll have to dig a bit deeper  to
 see what's going on.

 Dmitry.


 On Mon Feb 23 2015 at 10:03:33 AM David Knezevic 
 david.kneze...@akselos.com wrote:

 Hi Dmitry,

 OK, good to hear we're seeing the same behavior for the example.

 Regarding this comment:


 libMesh needs to change the way configure extracts PETSc information --
 configuration data were moved:
 conf -- lib/petsc-conf
 ${PETSC_ARCH}/conf -- ${PETSC_ARCH}/lib/petsc-conf

 At one point I started looking at m4/petsc.m4, but that got put on the
 back burner.  For now making the relevant symlinks by hand lets you
 configure and build libMesh with petsc/master.



 So do you suggest that the next step here is to build libmesh against
 petsc/master so that we can try your PETSc pull request that implements
 MatReset() to see if that gets this example working?

 David





 On Mon Feb 23 2015 at 9:15:44 AM David Knezevic 
 david.kneze...@akselos.com wrote:

 Hi Dmitry,

 Thanks very much for testing out the example.

 examples/systems_of_equations/ex8 works fine for me in serial, but
 it fails for me if I run with more than 1 MPI process. Can you try it 
 with,
 say, 2 or 4 MPI processes?

 If we need access to MatReset in libMesh to get this to work, I'll be
 happy to work on a libMesh pull request for that.

 David


 --

 David J. Knezevic | CTO
 Akselos | 17 Bay State Road | Boston, MA | 02215
 Phone (office): +1-857-265-2238
 Phone (mobile): +1-617-599-4755
 Web: http://www.akselos.com


 On Mon, Feb 23, 2015 at 10:08 AM, Dmitry Karpeyev 
 karp...@mcs.anl.gov wrote:

 David,

 What code are you running when you encounter this error?  I'm trying
 to reproduce it and
 I tried examples/systems_of_equations/ex8, but it ran for me,
 ostensibly to completion.

 I have a small PETSc pull request that implements MatReset(), which
 passes a small PETSc test,
 but libMesh needs some work to be able to build against petsc/master
 because of some recent
 changes to PETSc.

 Dmitry.

 On Mon Feb 23 2015 at 7:17:06 AM David Knezevic 
 david.kneze...@akselos.com wrote:

 Hi Barry, hi Dmitry,

 I set the matrix to BAIJ and back to AIJ, and the code got a bit
 further. But I now run into the error pasted below (Note that I'm now 
 using
 --with-debugging=1):

 PETSC ERROR: - Error Message
 --
 PETSC ERROR: Petsc has generated inconsistent data
 PETSC ERROR: MPIAIJ Matrix was assembled but is missing

Re: [petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length

[David Knezevic]

Hi Barry, hi Dmitry,

I set the matrix to BAIJ and back to AIJ, and the code got a bit further.
But I now run into the error pasted below (Note that I'm now using
--with-debugging=1):

PETSC ERROR: - Error Message
--
PETSC ERROR: Petsc has generated inconsistent data
PETSC ERROR: MPIAIJ Matrix was assembled but is missing garray
PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014
PETSC ERROR: ./example-dbg on a arch-linux2-c-debug named david-Lenovo by
dknez Mon Feb 23 08:05:44 2015
PETSC ERROR: Configure options --with-shared-libraries=1 --with-debugging=1
--download-suitesparse=1 --download-parmetis=1 --download-blacs=1
--download-scalapack=1 --download-mumps=1 --download-metis
--download-superlu_dist
--prefix=/home/dknez/software/libmesh_install/dbg_real/petsc
--download-hypre
PETSC ERROR: #1 MatCreateColmap_MPIAIJ_Private() line 361 in
/home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
PETSC ERROR: #2 MatSetValues_MPIAIJ() line 538 in
/home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
PETSC ERROR: #3 MatSetValues() line 1136 in
/home/dknez/software/petsc-3.5.2/src/mat/interface/matrix.c
PETSC ERROR: #4 add_matrix() line 765 in
/home/dknez/software/libmesh-src/src/numerics/petsc_matrix.C
--

This occurs when I try to set some entries of the matrix. Do you have any
suggestions on how I can resolve this?

Thanks!
David




On Sun, Feb 22, 2015 at 10:22 PM, Dmitry Karpeyev dkarp...@gmail.com
wrote:



 On Sun Feb 22 2015 at 9:15:22 PM Barry Smith bsm...@mcs.anl.gov wrote:


  On Feb 22, 2015, at 9:09 PM, David Knezevic david.kneze...@akselos.com
 wrote:
 
  Hi Dmitry,
 
  Thanks for the suggestion. I tried MatSetType(mat,MATMPIAIJ) followed
 by MatXAIJSetPreallocation(...), but unfortunately this still gives me the
 same error as before: nnz cannot be greater than row length: local row 168
 value 24 rowlength 0.
 
  I gather that the idea here is that MatSetType builds a new matrix
 object, and then I should be able to pre-allocate for that new matrix
 however I like, right? Was I supposed to clear the matrix object somehow
 before calling MatSetType? (I didn't do any sort of clear operation.)

   If the type doesn't change then MatSetType() won't do anything. You can
 try setting the type to BAIJ and then setting the type back to AIJ. This
 may/should clear out the matrix.

 Ah, yes.  If the type is the same as before it does quit early, but
 changing the type and then back will clear out and rebuild the matrix.   We
 need
 something like MatReset() to do the equivalent thing.


 
  As I said earlier, I'll make a dbg PETSc build, so hopefully that will
 help shed some light on what's going wrong for me.

 I think it's always a good idea to have a dbg build of PETSc when you
 doing things like these.

 Dmitry.


Don't bother, what I suggested won't work.

   Barry


 
  Thanks,
  David
 
 
 
 
  On Sun, Feb 22, 2015 at 6:02 PM, Dmitry Karpeyev dkarp...@gmail.com
 wrote:
  David,
  It might be easier to just rebuild the whole matrix from scratch: you
 would in effect be doing all that with disassembling and resetting the
 preallocation.
  MatSetType(mat,MATMPIAIJ)
  or
  PetscObjectGetType((PetscObject)mat,type);
  MatSetType(mat,type);
  followed by
  MatXAIJSetPreallocation(...);
  should do.
  Dmitry.
 
 
  On Sun Feb 22 2015 at 4:45:46 PM Barry Smith bsm...@mcs.anl.gov
 wrote:
 
   Do not call for SeqAIJ matrix. Do not call before the first time you
 have preallocated and put entries in the matrix and done the
 MatAssemblyBegin/End()
 
If it still crashes you'll need to try the debugger
 
Barry
 
   On Feb 22, 2015, at 4:09 PM, David Knezevic 
 david.kneze...@akselos.com wrote:
  
   Hi Barry,
  
   Thanks for your help, much appreciated.
  
   I added a prototype for MatDisAssemble_MPIAIJ:
   PETSC_INTERN PetscErrorCode MatDisAssemble_MPIAIJ(Mat);
  
   and I added a call to MatDisAssemble_MPIAIJ before
 MatMPIAIJSetPreallocation. However, I get a segfault on the call to
 MatDisAssemble_MPIAIJ. The segfault occurs in both serial and parallel.
  
   FYI, I'm using Petsc 3.5.2, and I'm not using a non-debug build
 (though I could rebuild PETSc in debug mode if you think that would help
 figure out what's happening here).
  
   Thanks,
   David
  
  
  
   On Sun, Feb 22, 2015 at 1:13 PM, Barry Smith bsm...@mcs.anl.gov
 wrote:
  
 David,
  
  This is an obscure little feature of MatMPIAIJ,   each time you
 change the sparsity pattern before you call the MatMPIAIJSetPreallocation
 you need to call  MatDisAssemble_MPIAIJ(Mat mat).This is a private
 PETSc function so you need to provide your own prototype for it above the
 function you use it in.
  
 Let us know if this resolves the problem.
  
  Barry

Re: [petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length

[David Knezevic]

Hi Dmitry,

Thanks very much for testing out the example.

examples/systems_of_equations/ex8 works fine for me in serial, but it fails
for me if I run with more than 1 MPI process. Can you try it with, say, 2
or 4 MPI processes?

If we need access to MatReset in libMesh to get this to work, I'll be happy
to work on a libMesh pull request for that.

David


--

David J. Knezevic | CTO
Akselos | 17 Bay State Road | Boston, MA | 02215
Phone (office): +1-857-265-2238
Phone (mobile): +1-617-599-4755
Web: http://www.akselos.com


On Mon, Feb 23, 2015 at 10:08 AM, Dmitry Karpeyev karp...@mcs.anl.gov
wrote:

 David,

 What code are you running when you encounter this error?  I'm trying to
 reproduce it and
 I tried examples/systems_of_equations/ex8, but it ran for me, ostensibly
 to completion.

 I have a small PETSc pull request that implements MatReset(), which passes
 a small PETSc test,
 but libMesh needs some work to be able to build against petsc/master
 because of some recent
 changes to PETSc.

 Dmitry.

 On Mon Feb 23 2015 at 7:17:06 AM David Knezevic 
 david.kneze...@akselos.com wrote:

 Hi Barry, hi Dmitry,

 I set the matrix to BAIJ and back to AIJ, and the code got a bit further.
 But I now run into the error pasted below (Note that I'm now using
 --with-debugging=1):

 PETSC ERROR: - Error Message
 --
 PETSC ERROR: Petsc has generated inconsistent data
 PETSC ERROR: MPIAIJ Matrix was assembled but is missing garray
 PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
 trouble shooting.
 PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014
 PETSC ERROR: ./example-dbg on a arch-linux2-c-debug named david-Lenovo by
 dknez Mon Feb 23 08:05:44 2015
 PETSC ERROR: Configure options --with-shared-libraries=1
 --with-debugging=1 --download-suitesparse=1 --download-parmetis=1
 --download-blacs=1 --download-scalapack=1 --download-mumps=1
 --download-metis --download-superlu_dist
 --prefix=/home/dknez/software/libmesh_install/dbg_real/petsc
 --download-hypre
 PETSC ERROR: #1 MatCreateColmap_MPIAIJ_Private() line 361 in
 /home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
 PETSC ERROR: #2 MatSetValues_MPIAIJ() line 538 in
 /home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
 PETSC ERROR: #3 MatSetValues() line 1136 in
 /home/dknez/software/petsc-3.5.2/src/mat/interface/matrix.c
 PETSC ERROR: #4 add_matrix() line 765 in
 /home/dknez/software/libmesh-src/src/numerics/petsc_matrix.C
 --

 This occurs when I try to set some entries of the matrix. Do you have any
 suggestions on how I can resolve this?

 Thanks!
 David




 On Sun, Feb 22, 2015 at 10:22 PM, Dmitry Karpeyev dkarp...@gmail.com
 wrote:



 On Sun Feb 22 2015 at 9:15:22 PM Barry Smith bsm...@mcs.anl.gov wrote:


  On Feb 22, 2015, at 9:09 PM, David Knezevic 
 david.kneze...@akselos.com wrote:
 
  Hi Dmitry,
 
  Thanks for the suggestion. I tried MatSetType(mat,MATMPIAIJ) followed
 by MatXAIJSetPreallocation(...), but unfortunately this still gives me the
 same error as before: nnz cannot be greater than row length: local row 168
 value 24 rowlength 0.
 
  I gather that the idea here is that MatSetType builds a new matrix
 object, and then I should be able to pre-allocate for that new matrix
 however I like, right? Was I supposed to clear the matrix object somehow
 before calling MatSetType? (I didn't do any sort of clear operation.)

   If the type doesn't change then MatSetType() won't do anything. You
 can try setting the type to BAIJ and then setting the type back to AIJ.
 This may/should clear out the matrix.

 Ah, yes.  If the type is the same as before it does quit early, but
 changing the type and then back will clear out and rebuild the matrix.   We
 need
 something like MatReset() to do the equivalent thing.


 
  As I said earlier, I'll make a dbg PETSc build, so hopefully that
 will help shed some light on what's going wrong for me.

 I think it's always a good idea to have a dbg build of PETSc when you
 doing things like these.

 Dmitry.


Don't bother, what I suggested won't work.

   Barry


 
  Thanks,
  David
 
 
 
 
  On Sun, Feb 22, 2015 at 6:02 PM, Dmitry Karpeyev dkarp...@gmail.com
 wrote:
  David,
  It might be easier to just rebuild the whole matrix from scratch: you
 would in effect be doing all that with disassembling and resetting the
 preallocation.
  MatSetType(mat,MATMPIAIJ)
  or
  PetscObjectGetType((PetscObject)mat,type);
  MatSetType(mat,type);
  followed by
  MatXAIJSetPreallocation(...);
  should do.
  Dmitry.
 
 
  On Sun Feb 22 2015 at 4:45:46 PM Barry Smith bsm...@mcs.anl.gov
 wrote:
 
   Do not call for SeqAIJ matrix. Do not call before the first time you
 have preallocated and put entries in the matrix and done the
 MatAssemblyBegin/End()
 
If it still crashes you'll need to try the debugger
 
Barry

Re: [petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length

[David Knezevic]

OK, sounds good. Let me know if I can help with the digging.

David



On Mon, Feb 23, 2015 at 11:10 AM, Dmitry Karpeyev karp...@mcs.anl.gov
wrote:

 I just tried building against petsc/master, but there needs to be more
 work on libMesh before it can work with petsc/master:
 the new VecLockPush()/Pop() stuff isn't respected by vector manipulation
 in libMesh.
 I put a hack equivalent to MatReset() into your branch (patch attached, in
 case you want to take a look at it),
 but it generates the same error in MatCreateColmap that you reported
 earlier.  It's odd that it occurs
 on the second nonlinear iteration.   I'll have to dig a bit deeper  to
 see what's going on.

 Dmitry.


 On Mon Feb 23 2015 at 10:03:33 AM David Knezevic 
 david.kneze...@akselos.com wrote:

 Hi Dmitry,

 OK, good to hear we're seeing the same behavior for the example.

 Regarding this comment:


 libMesh needs to change the way configure extracts PETSc information --
 configuration data were moved:
 conf -- lib/petsc-conf
 ${PETSC_ARCH}/conf -- ${PETSC_ARCH}/lib/petsc-conf

 At one point I started looking at m4/petsc.m4, but that got put on the
 back burner.  For now making the relevant symlinks by hand lets you
 configure and build libMesh with petsc/master.



 So do you suggest that the next step here is to build libmesh against
 petsc/master so that we can try your PETSc pull request that implements
 MatReset() to see if that gets this example working?

 David





 On Mon Feb 23 2015 at 9:15:44 AM David Knezevic 
 david.kneze...@akselos.com wrote:

 Hi Dmitry,

 Thanks very much for testing out the example.

 examples/systems_of_equations/ex8 works fine for me in serial, but it
 fails for me if I run with more than 1 MPI process. Can you try it with,
 say, 2 or 4 MPI processes?

 If we need access to MatReset in libMesh to get this to work, I'll be
 happy to work on a libMesh pull request for that.

 David


 --

 David J. Knezevic | CTO
 Akselos | 17 Bay State Road | Boston, MA | 02215
 Phone (office): +1-857-265-2238
 Phone (mobile): +1-617-599-4755
 Web: http://www.akselos.com


 On Mon, Feb 23, 2015 at 10:08 AM, Dmitry Karpeyev karp...@mcs.anl.gov
 wrote:

 David,

 What code are you running when you encounter this error?  I'm trying
 to reproduce it and
 I tried examples/systems_of_equations/ex8, but it ran for me,
 ostensibly to completion.

 I have a small PETSc pull request that implements MatReset(), which
 passes a small PETSc test,
 but libMesh needs some work to be able to build against petsc/master
 because of some recent
 changes to PETSc.

 Dmitry.

 On Mon Feb 23 2015 at 7:17:06 AM David Knezevic 
 david.kneze...@akselos.com wrote:

 Hi Barry, hi Dmitry,

 I set the matrix to BAIJ and back to AIJ, and the code got a bit
 further. But I now run into the error pasted below (Note that I'm now 
 using
 --with-debugging=1):

 PETSC ERROR: - Error Message
 --
 PETSC ERROR: Petsc has generated inconsistent data
 PETSC ERROR: MPIAIJ Matrix was assembled but is missing garray
 PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
 for trouble shooting.
 PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014
 PETSC ERROR: ./example-dbg on a arch-linux2-c-debug named
 david-Lenovo by dknez Mon Feb 23 08:05:44 2015
 PETSC ERROR: Configure options --with-shared-libraries=1
 --with-debugging=1 --download-suitesparse=1 --download-parmetis=1
 --download-blacs=1 --download-scalapack=1 --download-mumps=1
 --download-metis --download-superlu_dist 
 --prefix=/home/dknez/software/libmesh_install/dbg_real/petsc
 --download-hypre
 PETSC ERROR: #1 MatCreateColmap_MPIAIJ_Private() line 361 in
 /home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
 PETSC ERROR: #2 MatSetValues_MPIAIJ() line 538 in
 /home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
 PETSC ERROR: #3 MatSetValues() line 1136 in
 /home/dknez/software/petsc-3.5.2/src/mat/interface/matrix.c
 PETSC ERROR: #4 add_matrix() line 765 in /home/dknez/software/libmesh-
 src/src/numerics/petsc_matrix.C
 
 --

 This occurs when I try to set some entries of the matrix. Do you have
 any suggestions on how I can resolve this?

 Thanks!
 David




 On Sun, Feb 22, 2015 at 10:22 PM, Dmitry Karpeyev dkarp...@gmail.com
  wrote:



 On Sun Feb 22 2015 at 9:15:22 PM Barry Smith bsm...@mcs.anl.gov
 wrote:


  On Feb 22, 2015, at 9:09 PM, David Knezevic 
 david.kneze...@akselos.com wrote:
 
  Hi Dmitry,
 
  Thanks for the suggestion. I tried MatSetType(mat,MATMPIAIJ)
 followed by MatXAIJSetPreallocation(...), but unfortunately this still
 gives me the same error as before: nnz cannot be greater than row 
 length:
 local row 168 value 24 rowlength 0.
 
  I gather that the idea here is that MatSetType builds a new
 matrix object, and then I should be able to pre-allocate for that new
 matrix however I like

Re: [petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length

[David Knezevic]

Hi Dmitry,

OK, good to hear we're seeing the same behavior for the example.

Regarding this comment:


libMesh needs to change the way configure extracts PETSc information --
 configuration data were moved:
 conf -- lib/petsc-conf
 ${PETSC_ARCH}/conf -- ${PETSC_ARCH}/lib/petsc-conf

 At one point I started looking at m4/petsc.m4, but that got put on the
 back burner.  For now making the relevant symlinks by hand lets you
 configure and build libMesh with petsc/master.



So do you suggest that the next step here is to build libmesh against
petsc/master so that we can try your PETSc pull request that implements
MatReset() to see if that gets this example working?

David





 On Mon Feb 23 2015 at 9:15:44 AM David Knezevic 
 david.kneze...@akselos.com wrote:

 Hi Dmitry,

 Thanks very much for testing out the example.

 examples/systems_of_equations/ex8 works fine for me in serial, but it
 fails for me if I run with more than 1 MPI process. Can you try it with,
 say, 2 or 4 MPI processes?

 If we need access to MatReset in libMesh to get this to work, I'll be
 happy to work on a libMesh pull request for that.

 David


 --

 David J. Knezevic | CTO
 Akselos | 17 Bay State Road | Boston, MA | 02215
 Phone (office): +1-857-265-2238
 Phone (mobile): +1-617-599-4755
 Web: http://www.akselos.com


 On Mon, Feb 23, 2015 at 10:08 AM, Dmitry Karpeyev karp...@mcs.anl.gov
 wrote:

 David,

 What code are you running when you encounter this error?  I'm trying to
 reproduce it and
 I tried examples/systems_of_equations/ex8, but it ran for me,
 ostensibly to completion.

 I have a small PETSc pull request that implements MatReset(), which
 passes a small PETSc test,
 but libMesh needs some work to be able to build against petsc/master
 because of some recent
 changes to PETSc.

 Dmitry.

 On Mon Feb 23 2015 at 7:17:06 AM David Knezevic 
 david.kneze...@akselos.com wrote:

 Hi Barry, hi Dmitry,

 I set the matrix to BAIJ and back to AIJ, and the code got a bit
 further. But I now run into the error pasted below (Note that I'm now using
 --with-debugging=1):

 PETSC ERROR: - Error Message
 --
 PETSC ERROR: Petsc has generated inconsistent data
 PETSC ERROR: MPIAIJ Matrix was assembled but is missing garray
 PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
 for trouble shooting.
 PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014
 PETSC ERROR: ./example-dbg on a arch-linux2-c-debug named david-Lenovo
 by dknez Mon Feb 23 08:05:44 2015
 PETSC ERROR: Configure options --with-shared-libraries=1
 --with-debugging=1 --download-suitesparse=1 --download-parmetis=1
 --download-blacs=1 --download-scalapack=1 --download-mumps=1
 --download-metis --download-superlu_dist 
 --prefix=/home/dknez/software/libmesh_install/dbg_real/petsc
 --download-hypre
 PETSC ERROR: #1 MatCreateColmap_MPIAIJ_Private() line 361 in
 /home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
 PETSC ERROR: #2 MatSetValues_MPIAIJ() line 538 in
 /home/dknez/software/petsc-3.5.2/src/mat/impls/aij/mpi/mpiaij.c
 PETSC ERROR: #3 MatSetValues() line 1136 in
 /home/dknez/software/petsc-3.5.2/src/mat/interface/matrix.c
 PETSC ERROR: #4 add_matrix() line 765 in /home/dknez/software/libmesh-
 src/src/numerics/petsc_matrix.C
 
 --

 This occurs when I try to set some entries of the matrix. Do you have
 any suggestions on how I can resolve this?

 Thanks!
 David




 On Sun, Feb 22, 2015 at 10:22 PM, Dmitry Karpeyev dkarp...@gmail.com
 wrote:



 On Sun Feb 22 2015 at 9:15:22 PM Barry Smith bsm...@mcs.anl.gov
 wrote:


  On Feb 22, 2015, at 9:09 PM, David Knezevic 
 david.kneze...@akselos.com wrote:
 
  Hi Dmitry,
 
  Thanks for the suggestion. I tried MatSetType(mat,MATMPIAIJ)
 followed by MatXAIJSetPreallocation(...), but unfortunately this still
 gives me the same error as before: nnz cannot be greater than row 
 length:
 local row 168 value 24 rowlength 0.
 
  I gather that the idea here is that MatSetType builds a new matrix
 object, and then I should be able to pre-allocate for that new matrix
 however I like, right? Was I supposed to clear the matrix object somehow
 before calling MatSetType? (I didn't do any sort of clear operation.)

   If the type doesn't change then MatSetType() won't do anything. You
 can try setting the type to BAIJ and then setting the type back to AIJ.
 This may/should clear out the matrix.

 Ah, yes.  If the type is the same as before it does quit early, but
 changing the type and then back will clear out and rebuild the matrix.   
 We
 need
 something like MatReset() to do the equivalent thing.


 
  As I said earlier, I'll make a dbg PETSc build, so hopefully that
 will help shed some light on what's going wrong for me.

 I think it's always a good idea to have a dbg build of PETSc when you
 doing things like these.

 Dmitry.


Don't bother, what I

[petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length

[David Knezevic]

Hi all,

I've implemented a solver for a contact problem using SNES. The sparsity
pattern of the jacobian matrix needs to change at each nonlinear iteration
(because the elements which are in contact can change), so I tried to deal
with this by calling MatSeqAIJSetPreallocation
and MatMPIAIJSetPreallocation during each iteration in order to update the
preallocation.

This seems to work fine in serial, but with two or more MPI processes I run
into the error nnz cannot be greater than row length, e.g.:
nnz cannot be greater than row length: local row 528 value 12 rowlength 0

This error is from the call to
MatSeqAIJSetPreallocation(b-B,o_nz,o_nnz);
in MatMPIAIJSetPreallocation_MPIAIJ.

Any guidance on what the problem might be would be most appreciated. For
example, I was wondering if there is a problem with calling
SetPreallocation on a matrix that has already been preallocated?

Some notes:
- I'm using PETSc via libMesh
- The code that triggers this issue is available as a PR on the libMesh
github repo, in case anyone is interested:
https://github.com/libMesh/libmesh/pull/460/
- I can try to make a minimal pure-PETSc example that reproduces this
error, if that would be helpful.

Many thanks,
David

Re: [petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length

[David Knezevic]

Thanks, that helps! After fixing that, now I get this error:

[1]PETSC ERROR: Petsc has generated inconsistent data
[1]PETSC ERROR: MPIAIJ Matrix was assembled but is missing garray

Any suggestions about what may be wrong now? I'll try the debugger tomorrow.

Thanks,
David



On Sun, Feb 22, 2015 at 5:45 PM, Barry Smith bsm...@mcs.anl.gov wrote:


  Do not call for SeqAIJ matrix. Do not call before the first time you have
 preallocated and put entries in the matrix and done the
 MatAssemblyBegin/End()

   If it still crashes you'll need to try the debugger

   Barry

  On Feb 22, 2015, at 4:09 PM, David Knezevic david.kneze...@akselos.com
 wrote:
 
  Hi Barry,
 
  Thanks for your help, much appreciated.
 
  I added a prototype for MatDisAssemble_MPIAIJ:
  PETSC_INTERN PetscErrorCode MatDisAssemble_MPIAIJ(Mat);
 
  and I added a call to MatDisAssemble_MPIAIJ before
 MatMPIAIJSetPreallocation. However, I get a segfault on the call to
 MatDisAssemble_MPIAIJ. The segfault occurs in both serial and parallel.
 
  FYI, I'm using Petsc 3.5.2, and I'm not using a non-debug build (though
 I could rebuild PETSc in debug mode if you think that would help figure out
 what's happening here).
 
  Thanks,
  David
 
 
 
  On Sun, Feb 22, 2015 at 1:13 PM, Barry Smith bsm...@mcs.anl.gov wrote:
 
David,
 
 This is an obscure little feature of MatMPIAIJ,   each time you
 change the sparsity pattern before you call the MatMPIAIJSetPreallocation
 you need to call  MatDisAssemble_MPIAIJ(Mat mat).This is a private
 PETSc function so you need to provide your own prototype for it above the
 function you use it in.
 
Let us know if this resolves the problem.
 
 Barry
 
  We never really intended that people would call
 MatMPIAIJSetPreallocation() AFTER they had already used the matrix.
 
 
   On Feb 22, 2015, at 6:50 AM, David Knezevic 
 david.kneze...@akselos.com wrote:
  
   Hi all,
  
   I've implemented a solver for a contact problem using SNES. The
 sparsity pattern of the jacobian matrix needs to change at each nonlinear
 iteration (because the elements which are in contact can change), so I
 tried to deal with this by calling MatSeqAIJSetPreallocation and
 MatMPIAIJSetPreallocation during each iteration in order to update the
 preallocation.
  
   This seems to work fine in serial, but with two or more MPI processes
 I run into the error nnz cannot be greater than row length, e.g.:
   nnz cannot be greater than row length: local row 528 value 12
 rowlength 0
  
   This error is from the call to
   MatSeqAIJSetPreallocation(b-B,o_nz,o_nnz); in
 MatMPIAIJSetPreallocation_MPIAIJ.
  
   Any guidance on what the problem might be would be most appreciated.
 For example, I was wondering if there is a problem with calling
 SetPreallocation on a matrix that has already been preallocated?
  
   Some notes:
   - I'm using PETSc via libMesh
   - The code that triggers this issue is available as a PR on the
 libMesh github repo, in case anyone is interested:
 https://github.com/libMesh/libmesh/pull/460/
   - I can try to make a minimal pure-PETSc example that reproduces this
 error, if that would be helpful.
  
   Many thanks,
   David
  
 
 

Re: [petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length

[David Knezevic]

P.S. Typo correction: I meant to say I'm using a non-debug build.

David



On Sun, Feb 22, 2015 at 5:09 PM, David Knezevic david.kneze...@akselos.com
wrote:

 Hi Barry,

 Thanks for your help, much appreciated.

 I added a prototype for MatDisAssemble_MPIAIJ:
 PETSC_INTERN PetscErrorCode MatDisAssemble_MPIAIJ(Mat);

 and I added a call to MatDisAssemble_MPIAIJ before
 MatMPIAIJSetPreallocation. However, I get a segfault on the call
 to MatDisAssemble_MPIAIJ. The segfault occurs in both serial and parallel.

 FYI, I'm using Petsc 3.5.2, and I'm not using a non-debug build (though I
 could rebuild PETSc in debug mode if you think that would help figure out
 what's happening here).

 Thanks,
 David



 On Sun, Feb 22, 2015 at 1:13 PM, Barry Smith bsm...@mcs.anl.gov wrote:


   David,

This is an obscure little feature of MatMPIAIJ,   each time you change
 the sparsity pattern before you call the MatMPIAIJSetPreallocation you need
 to call  MatDisAssemble_MPIAIJ(Mat mat).This is a private PETSc
 function so you need to provide your own prototype for it above the
 function you use it in.

   Let us know if this resolves the problem.

Barry

 We never really intended that people would call
 MatMPIAIJSetPreallocation() AFTER they had already used the matrix.


  On Feb 22, 2015, at 6:50 AM, David Knezevic david.kneze...@akselos.com
 wrote:
 
  Hi all,
 
  I've implemented a solver for a contact problem using SNES. The
 sparsity pattern of the jacobian matrix needs to change at each nonlinear
 iteration (because the elements which are in contact can change), so I
 tried to deal with this by calling MatSeqAIJSetPreallocation and
 MatMPIAIJSetPreallocation during each iteration in order to update the
 preallocation.
 
  This seems to work fine in serial, but with two or more MPI processes I
 run into the error nnz cannot be greater than row length, e.g.:
  nnz cannot be greater than row length: local row 528 value 12 rowlength
 0
 
  This error is from the call to
  MatSeqAIJSetPreallocation(b-B,o_nz,o_nnz); in
 MatMPIAIJSetPreallocation_MPIAIJ.
 
  Any guidance on what the problem might be would be most appreciated.
 For example, I was wondering if there is a problem with calling
 SetPreallocation on a matrix that has already been preallocated?
 
  Some notes:
  - I'm using PETSc via libMesh
  - The code that triggers this issue is available as a PR on the libMesh
 github repo, in case anyone is interested:
 https://github.com/libMesh/libmesh/pull/460/
  - I can try to make a minimal pure-PETSc example that reproduces this
 error, if that would be helpful.
 
  Many thanks,
  David
 

Re: [petsc-users] MatMPIAIJSetPreallocation: nnz cannot be greater than row length

[David Knezevic]

Hi Dmitry,

Thanks for the suggestion. I tried MatSetType(mat,MATMPIAIJ) followed
by MatXAIJSetPreallocation(...),
but unfortunately this still gives me the same error as before: nnz cannot
be greater than row length: local row 168 value 24 rowlength 0.

I gather that the idea here is that MatSetType builds a new matrix object,
and then I should be able to pre-allocate for that new matrix however I
like, right? Was I supposed to clear the matrix object somehow before
calling MatSetType? (I didn't do any sort of clear operation.)

As I said earlier, I'll make a dbg PETSc build, so hopefully that will help
shed some light on what's going wrong for me.

Thanks,
David




On Sun, Feb 22, 2015 at 6:02 PM, Dmitry Karpeyev dkarp...@gmail.com wrote:

 David,
 It might be easier to just rebuild the whole matrix from scratch: you
 would in effect be doing all that with disassembling and resetting the
 preallocation.
 MatSetType(mat,MATMPIAIJ)
 or
 PetscObjectGetType((PetscObject)mat,type);
 MatSetType(mat,type);
 followed by
 MatXAIJSetPreallocation(...);
 should do.
 Dmitry.


 On Sun Feb 22 2015 at 4:45:46 PM Barry Smith bsm...@mcs.anl.gov wrote:


  Do not call for SeqAIJ matrix. Do not call before the first time you
 have preallocated and put entries in the matrix and done the
 MatAssemblyBegin/End()

   If it still crashes you'll need to try the debugger

   Barry

  On Feb 22, 2015, at 4:09 PM, David Knezevic david.kneze...@akselos.com
 wrote:
 
  Hi Barry,
 
  Thanks for your help, much appreciated.
 
  I added a prototype for MatDisAssemble_MPIAIJ:
  PETSC_INTERN PetscErrorCode MatDisAssemble_MPIAIJ(Mat);
 
  and I added a call to MatDisAssemble_MPIAIJ before
 MatMPIAIJSetPreallocation. However, I get a segfault on the call to
 MatDisAssemble_MPIAIJ. The segfault occurs in both serial and parallel.
 
  FYI, I'm using Petsc 3.5.2, and I'm not using a non-debug build (though
 I could rebuild PETSc in debug mode if you think that would help figure out
 what's happening here).
 
  Thanks,
  David
 
 
 
  On Sun, Feb 22, 2015 at 1:13 PM, Barry Smith bsm...@mcs.anl.gov
 wrote:
 
David,
 
 This is an obscure little feature of MatMPIAIJ,   each time you
 change the sparsity pattern before you call the MatMPIAIJSetPreallocation
 you need to call  MatDisAssemble_MPIAIJ(Mat mat).This is a private
 PETSc function so you need to provide your own prototype for it above the
 function you use it in.
 
Let us know if this resolves the problem.
 
 Barry
 
  We never really intended that people would call
 MatMPIAIJSetPreallocation() AFTER they had already used the matrix.
 
 
   On Feb 22, 2015, at 6:50 AM, David Knezevic 
 david.kneze...@akselos.com wrote:
  
   Hi all,
  
   I've implemented a solver for a contact problem using SNES. The
 sparsity pattern of the jacobian matrix needs to change at each nonlinear
 iteration (because the elements which are in contact can change), so I
 tried to deal with this by calling MatSeqAIJSetPreallocation and
 MatMPIAIJSetPreallocation during each iteration in order to update the
 preallocation.
  
   This seems to work fine in serial, but with two or more MPI processes
 I run into the error nnz cannot be greater than row length, e.g.:
   nnz cannot be greater than row length: local row 528 value 12
 rowlength 0
  
   This error is from the call to
   MatSeqAIJSetPreallocation(b-B,o_nz,o_nnz); in
 MatMPIAIJSetPreallocation_MPIAIJ.
  
   Any guidance on what the problem might be would be most appreciated.
 For example, I was wondering if there is a problem with calling
 SetPreallocation on a matrix that has already been preallocated?
  
   Some notes:
   - I'm using PETSc via libMesh
   - The code that triggers this issue is available as a PR on the
 libMesh github repo, in case anyone is interested:
 https://github.com/libMesh/libmesh/pull/460/
   - I can try to make a minimal pure-PETSc example that reproduces this
 error, if that would be helpful.
  
   Many thanks,
   David
  
 
 

Stalling once linear system becomes a certain size

[David Knezevic]

Hello,

I am trying to run a PETSc code on a parallel machine (it may be 
relevant that each node contains four AMD Opteron Quad-Core 64-bit 
processors (16 cores in all) as an SMP unit with 32GB of memory) and I'm 
observing some behaviour I don't understand.

I'm using PETSC_COMM_SELF in order to construct the same matrix on each 
processor (and solve the system with a different right-hand side vector 
on each processor), and when each linear system is around 315x315 
(block-sparse), then each linear system is solved very quickly on each 
processor (approx 7x10^{-4} seconds), but when I increase the size of 
the linear system to 350x350 (or larger), the linear solves completely 
stall. I've tried a number of different solvers and preconditioners, but 
nothing seems to help. Also, this code has worked very well on other 
machines, although the machines I have used it on before have not had 
this architecture in which each node is an SMP unit. I was wondering if 
you have observed this kind of issue before?

I'm using PETSc 2.3.3, compiled with the Intel 10.1 compiler.

Thanks very much,
David

Using SUBSET_NONZERO_PATTERN with MatAXPY

[David Knezevic]

Hello,

I've got a SeqBAIJ matrix A which I am reassembling at each time-step in 
order to solve a time-dependent PDE. I initially set up the sparsity 
pattern of A, and then assemble A at each time step by doing operations 
like the following:

ierr = MatZeroEntries(A);CHKERRQ(ierr);
ierr = MatAXPY(A,const1,B1,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = MatAXPY(A,const2,B2,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);

This works fine. However, the sparsity patterns of B1, B2 etc should all 
be subsets of the sparsity pattern of A (which should be retained by 
MatZeroEntries). But changing DIFFERENT_NONZERO_PATTERN to 
SUBSET_NONZERO_PATTERN in the above code causes an error. In particular, 
the code runs fine for the first time step, but on the second time step 
I get the error:

[0]PETSC ERROR: - Error Message 

[0]PETSC ERROR: Invalid argument!
[0]PETSC ERROR: Wrong type of object: Parameter # 3!
...
[0]PETSC ERROR: MatAXPY() line 34 in src/mat/utils/axpy.c

I found a description of the same problem here:
http://ganymed.iwr.uni-heidelberg.de/pipermail/dealii/2007/002066.html
but I haven't found an answer. Assistance would be most appreciated.

Best regards,
David