[petsc-dev] coding style
Matthew Knepley
knepley at gmail.com
Thu Aug 18 09:00:54 CDT 2016
On Thu, Aug 18, 2016 at 8:55 AM, Munson, Todd <tmunson at mcs.anl.gov> wrote:
>
> Just to add some more thoughts to this argument. Currently TAO uses a
> model where we have methods to get the gradient and Hessian of the
> objective function.
>
> For quadratic problems we really want the Hessian matrix H and the
> vector c so that the objective function is x'*H*x + c'*x.
>
> Because we do not have any specialized code for quadratic problems,
> we end up using the gradient of the objective and the Hessian
> matrix to back out the constant c.
>
> c = grad - Hx
>
> where grad = Hx + c.
>
> This relies on cancellation of Hx in order to get an accurate value
> for the vector c, which can be a problem, especially if we start
> far from a solution.
>
What about c = grad 0?
Matt
> You may have the same problem with linear ODE if you are relying
> on a residual and jacobian evaluation pair to back out the
> constant right-hand side vector.
>
> Its probably too complicated for users, but what one might like is
> to have the user separate the objective function, for example, and
> provide, c, H, and a nonlinear part f. Then I can assemble the
> objective function internally as
>
> c'*x + x'*H*x + f(x)
>
> while retaining correct values for c and H when I need them.
>
> The user could continue as always by putting everything into f(x),
> but then they would have to use a nonlinear solver.
>
> I'd do the same for constraints:
>
> b + A*x + F(x)
>
> This also allows the methods to do internal type checking to ensure
> the problem passed into the methods has the correct signature.
> E.g. for quadratic programs the nonlinear functions are NULL.
>
> Todd.
>
> > On Aug 16, 2016, at 7:03 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> >
> >> On Aug 16, 2016, at 6:09 PM, Munson, Todd <tmunson at mcs.anl.gov> wrote:
> >>
> >>
> >> Is there a reason to put flags into PETSc for linear or nonlinear
> equations? In PETSc
> >> there is a natural difference between linear (KSP) and nonlinear (SNES)
> that is
> >> reflected in the objects.
> >
> > The flags are for TS ODE solvers (not KSP or SNES); they can come into
> play for two reasons
> >
> > 1) whether Jacobian entries can be changed inside the solver or not (for
> nonlinear ODES they can usually be changed since the user provides new
> values at each Jacobian update; for linear the user won't normally be
> putting in the same values again so if the solver changes the values it has
> to have some way to recover the values; often but not always the changes
> are only on the diagonal). Note that this reason is JUST a question of
> optimizing memory and code (do you make a copy of matrices or not?), not
> correctness of results.
> >
> > 2) There may also be specialized ODE integrator only worked for a linear
> problem. (I don't know if we have any of those implemented or not). This is
> not just a question of optimizing the solver, it is a question of
> correctness. If I try to solve a nonlinear problem with a linear only
> solver I will get a wrong or meaningless answer.
> >
> > Similarly I think for TAO both reasons for wanting to know properties
> of the problem exist.
> >>
> >> Any thoughts on splitting up the TAO methods in a similar way and
> distinguish at
> >> least between between quadratic and general nonlinear problems, rather
> than
> >> putting in flags. This would be consistent with the rest of PETSc.
> >>
> >> For example, instead of having a single "bound" directory, having
> something
> >> like a "bco" directory for generic nonlinear bound-constrained
> optimization
> >> and a "bqp" directory for bound constrained quadratic programs.
> >>
> >> I'd do the same for "unconstrained" into "uco" and "uqp".
> >>
> >> General constraints then become more cumbersome with too many
> combinations.
> >> Could probably get by with generic "lp", "qp", and "nlp" though.
> >
> > If you did this for the TAO user API then one could argue it should be
> done for TS also. But I don't think it should be done for TS.
> >
> > Two reasons for not having a separate API for linear and nonlinear
> problems is
> >
> > 1) I may want to use a "nonlinear" ODE integrator (that is an integrator
> that works for nonlinear problems also) on my linear problem and I
> shouldn't need to have ifs if my code to switch to use it. That is I don't
> want ugly code like
> >
> > if (using nonlinear integrator) {
> > TSSetIJacobian(..)
> > } else {
> > TSSetConstantMatrix(...)
> > }
> > Now I could also have my "nonlinear integrators" support the
> TSSetConstantMatrix() API but that leads to more complex implementations
> that are more work to maintain. Why not just have something like
> > TSSetIJacobian(ts, J,TSConstantJacobianFunction()) when you are
> doing a linear ODE?
> >
> > 2) I may spend 3 months writing my linear ODE problem using the linear
> interface and then realize I have to add a nonlinear term, now I have to
> change my app use the nonlinear interface of TS. This is why Matt argues we
> should just publicly support the SNES api even for linear problems.
> >
> > Dealing with changing or not changing matrices has another issue we
> have not dealt with properly in PETSc (or TAO). When the matrix entries are
> set. For linear ODE problems we kind of assume the matrix entries are set
> by the user before they ever call the TS code (and many users think this is
> a natural way to proceed), so essentially they compute the J matrix up
> front and then pass it to the TS which uses it. For nonlinear problems, of
> course, we cannot set the matrix entries up front, they MUST be set in the
> callback that TS calls to compute the matrix entries. But if the matrix
> structure is provided internally by a DM then the user does not have direct
> access to the matrix before calling TSSolve() so even for linear problems
> the code logic more naturally wants the (constant) Jacobian to be computed
> at the first timestep computation (and never again) rather than up front.
> Using the Matt Knepley consistency/simplicity argument we should just
> always use the callback approach and not even "support" the "linear only"
> approach. Then comes the nasty linear "flag" to indicate 1) if the user
> intends to provide the matrix only once and 2) to indicate the problem has
> special (linear) structure that only certain solvers will work correctly on.
> >
> > I'm rambling a bit but I need to in order to make clear these kind of
> design decisions are not straightforward.
> >
> > In terms of organizing the implementation code and the names (TaoType)
> of the methods I think spitting them into categories makes good sense. But
> never with your horrible abbreviations that no one but you would be able to
> remember! We are not using Fortran IV with 6 char limits.
> >
> > Before doing anything let's try to generate a table that includes all
> (most of?) the cases.
> >
> > The function being optimized is linear, quadratic, or general
> (another case is separable).
> >
> > The equality constraints are linear?, quadratic? general? or
> something else?
> >
> > The inequality constraints are simple bounds? linear? quadratic?
> general? something else?
> >
> > Is there a fourth (or even more) axis of problem properties?
> >
> > Is this level of categories enough to organize the directories and
> naming of TAO solvers?
> >
> > Barry
> >
> >
> >
> >>
> >> The methods would then explicitly mention the type of optimization
> problem
> >> with, for example, TAO_UCO_LMVM, TAO_BCO_TRON, and TAO_BQP_GPCG.
> >>
> >> Thoughts?
> >>
> >> Todd.
> >>
> >>
> >>> On Aug 16, 2016, at 4:32 PM, Oxberry, Geoffrey Malcolm <
> oxberry1 at llnl.gov> wrote:
> >>>
> >>> Definitely flags for linear problems would be helpful for TAO. Once
> there is an example up, I'd be happy to add that to the SQPTR pull request.
> >>>
> >>> ________________________________________
> >>> From: petsc-dev-bounces at mcs.anl.gov [petsc-dev-bounces at mcs.anl.gov]
> on behalf of Barry Smith [bsmith at mcs.anl.gov]
> >>> Sent: Monday, August 15, 2016 6:18 PM
> >>> To: Munson, Todd
> >>> Cc: petsc-dev
> >>> Subject: Re: [petsc-dev] coding style
> >>>
> >>>> On Aug 15, 2016, at 8:08 PM, Munson, Todd <tmunson at mcs.anl.gov>
> wrote:
> >>>>
> >>>>
> >>>> Since we are only doing diagonal modifications to the matrix in
> >>>> my case, should I simply create a matnest in my part of the
> >>>> code and apply the diagonal modification in the matnest
> >>>> without directly modifying any of the matrices?
> >>>> That might be cleaner for the users.
> >>>
> >>> It will be dependent on the solver how much of the matrix you need
> retain the parts you have modified.
> >>>
> >>> MatNest isn't the right thing here. I think you can do
> MatGetDiagonal() to keep the valid diagonal entries then MatDiagonalSet()
> to put back the valid entries.
> >>>
> >>> Barry
> >>>
> >>>
> >>>>
> >>>> Todd.
> >>>>
> >>>>> On Aug 15, 2016, at 7:55 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >>>>>
> >>>>>
> >>>>>> On Aug 15, 2016, at 7:42 PM, Munson, Todd <tmunson at mcs.anl.gov>
> wrote:
> >>>>>>
> >>>>>>
> >>>>>> Got it. I will get all the code fragments in the developer doc that
> >>>>>> are not meant to be deliberately wrong fixed tomorrow; at least
> >>>>>> that code will be compliant... :)
> >>>>>>
> >>>>>> Then its onto fixing bqpip.
> >>>>>>
> >>>>>> Was there a final say on the Hessian/Jacobian matrix when they are
> >>>>>> modified? Were we adding a "dirty" bit or was I changing my code
> >>>>>> to not modify those matrices?
> >>>>>
> >>>>> Changing email to petsc-dev since this is a general PETSc/TAO issue.
> >>>>>
> >>>>> This is why we don't wait months to fix something :-) I have
> completely forgot the context you are asking about.
> >>>>>
> >>>>> I think the general issue comes up when a problem is linear and the
> user just wants to set and forget the matrices while when the problem is
> nonlinear the user needs to reset values into the matrix anyways so it is
> fine for you to change the matrix internally. Unless there is a way for the
> user to explicitly indicate the system is linear your code cannot know if
> it needs to make a copy of the matrix.
> >>>>>
> >>>>> In TS we have typedef enum {TS_LINEAR,TS_NONLINEAR} TSProblemType;
> and TSSetProblemType() that can be used to indicate if the matrices can be
> changed or not (though I don't think we use it in this way, we have an open
> issue https://bitbucket.org/petsc/petsc/issues/135/memory-
> optimization-for-ts
> >>>>>
> >>>>> I think you should start by having a flag that indicates if it is
> safe to modify the matrices associated with the TAO object and if not safe
> keep a backup copy.
> >>>>>
> >>>>> Barry
> >>>>>
> >>>>>
> >>>>>
> >>>>>>
> >>>>>> Todd.
> >>>>>>
> >>>>>>> On Aug 15, 2016, at 6:09 PM, Barry Smith <bsmith at mcs.anl.gov>
> wrote:
> >>>>>>>
> >>>>>>>
> >>>>>>> Todd,
> >>>>>>>
> >>>>>>>> On Aug 15, 2016, at 10:26 AM, Satish Balay <balay at mcs.anl.gov>
> wrote:
> >>>>>>>>
> >>>>>>>> Todd,
> >>>>>>>>
> >>>>>>>> Forwarding this to Barry.
> >>>>>>>>
> >>>>>>>> Satish
> >>>>>>>>
> >>>>>>>> On Mon, 15 Aug 2016, Munson, Todd wrote:
> >>>>>>>>
> >>>>>>>>>
> >>>>>>>>> I'm going through the developer documentation. There seems to be
> >>>>>>>>> a few things missing from the coding style.
> >>>>>>>>>
> >>>>>>>>> For example, in function prototypes, is it:
> >>>>>>>>>
> >>>>>>>>> PetscErrorCode MyFunction(char *,const int *,int);
> >>>>>>>>>
> >>>>>>>>> or
> >>>>>>>>>
> >>>>>>>>> PetscErrorCode MyFunction(char*,const int*,int)
> >>>>>>>
> >>>>>>> The second one because the preference is to not have unnecessary
> space.
> >>>>>>>
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>> Same thing in the definitions:
> >>>>>>>>>
> >>>>>>>>> PetscErrorCode MyFunction(char *c,const int *i,int j)
> >>>>>>>>> {
> >>>>>>>>> }
> >>>>>>>>>
> >>>>>>>>> or
> >>>>>>>>>
> >>>>>>>>> PetscErrorCode MyFunction(char* c,const int* i,int j)
> >>>>>>>
> >>>>>>> The * should always be attached to the variable.
> >>>>>>>>> {
> >>>>>>>>> }
> >>>>>>>>>
> >>>>>>>>> The same questions apply with * replaced by [].
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>> Next is in the types for variables. Are pointers grouped
> >>>>>>>>> separately or together with the type:
> >>>>>>>>>
> >>>>>>>>> int *i,j,*k;
> >>>>>>>>>
> >>>>>>>>> or
> >>>>>>>>>
> >>>>>>>>> int *i,*k;
> >>>>>>>>> int j;
> >>>>>>>>>
> >>>>>>>>> If grouped together, should they be ordered with pointers first
> >>>>>>>>> and the others later?
> >>>>>>>>>
> >>>>>>>>> int *i,*k,j;
> >>>>>>>
> >>>>>>> There is no preference on the above; the pointers do not have to
> be on their own line.
> >>>>>>>>>
> >>>>>>>>> On the next topic, do you use
> >>>>>>>>>
> >>>>>>>>> char *blah;
> >>>>>>>>>
> >>>>>>>>> or
> >>>>>>>>>
> >>>>>>>>> char* blah;
> >>>>>>>>>
> >>>>>>>>> when there is only one of them?
> >>>>>>>
> >>>>>>> * belongs on the variable
> >>>>>>>
> >>>>>>>>>
> >>>>>>>>> I have the same issues with the typedefs. Sometimes there is
> >>>>>>>>>
> >>>>>>>>> typedef struct _EH* EH;
> >>>>>>>>>
> >>>>>>>>> and other times
> >>>>>>>>>
> >>>>>>>>> typedef struct _EH *EH;
> >>>>>>>
> >>>>>>> Here the * belongs as in the first case
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>> My personal thought is to be consistent and use char * with a
> space
> >>>>>>>>> everywhere, but you tell me. I will then go and fix the
> developer
> >>>>>>>>> documentation to be consistent with the style guide as a test of
> >>>>>>>>> my branch management.
> >>>>>>>>>
> >>>>>>>>> Todd.
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>
> >>>>>>>
> >>>>>>
> >>>>>
> >>>>
> >>>
> >>
> >
>
>
--
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20160818/3cddf32b/attachment.html>
More information about the petsc-dev
mailing list