> On Aug 16, 2016, at 6:09 PM, Munson, Todd <tmun...@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 <oxber...@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-boun...@mcs.anl.gov [petsc-dev-boun...@mcs.anl.gov] on
>> behalf of Barry Smith [bsm...@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 <tmun...@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 <bsm...@mcs.anl.gov> wrote:
>>>>
>>>>
>>>>> On Aug 15, 2016, at 7:42 PM, Munson, Todd <tmun...@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 <bsm...@mcs.anl.gov> wrote:
>>>>>>
>>>>>>
>>>>>> Todd,
>>>>>>
>>>>>>> On Aug 15, 2016, at 10:26 AM, Satish Balay <ba...@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.
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>
>>>>>
>>>>
>>>
>>
>
