Remove allocations for xl, xu, nlconstr0 in C interface

[nbelakovski]

As I was working on the Python bindings, I ran into some issues related to memory management. These changes should help with that, particularly the change to the allocation of nlconstr0, since that was the one that appeared to be giving issues with Python. In general I'd prefer to do as little memory allocation within the C interface as possible, since it's challenging enough to have memory management in 1 language, let alone 4 when we start interacting with Python->C++->C->Fortran.

What enables the removal of these allocations is that I've discovered that if a variable in Fortran is marked as a pointer and set to null(), then the Fortran function present returns false on that variable. This means we can use NULL from within C to indicate that a variable is not provided. This approach has a small inconvenience in that the user will have to handle allocations for optional variables.

For example, with f0, the user would have to do the following:

double f0 = obj_fun(x0);
problem->f0 = &f0;

Likewise for things like rhobeg, rhoend, the user would have to do something like the following:

double rhobeg = 1.0;
options->rhobeg = &rhobeg;
double rhoend = 1.0e-6;
options->rhoend = &rhoend;

This is in contrast to the current way of doing things:

options->rhobeg = 1.0;
options->rhoend = 1.0e-6;

As annoying as the previous option is, it might be preferable since it means we can stop duplicating initialization logic within C. That logic is currently implemented in prima_init_options and prima_check_problem, but it doesn't replicate the full logic for rhobeg/rhoend. I think that ideally C should be a 'dumb' pass-through and should let Fortran handle as much as possible.

We can use this PR to add this functionality to the rest of the variables, but I figured I would start with removing the memory allocations before going too much further.

[zaikunzhang]

Hi @nbelakovski , two quick questions:

1. Will this necessitate any change to the Fortran code?

2. For claims like the following, we need a reference to the standard.

   I've discovered that if a variable in Fortran is marked as a pointer and set to null(), then the Fortran function present returns false on that variable. This means we can use NULL from within C to indicate that a variable is not provided

Thank you.

[zaikunzhang]

It is the spirit of PRIMA to avoid requiring the user to do things like

double rhobeg = 1.0;
options->rhobeg = &rhobeg;
double rhoend = 1.0e-6;
options->rhoend = &rhoend;

This is unnatural and unfriendly. If this is inevitable, then we should add another layer to enable users to do options->rhobeg=1.0.

BTW, excuse me for my ignorance about C, could you remind me why it is an arrow instead of a dot?

N.B.: PRIMA tries to provide the most natural and friendly interface to users and handles everything else under the hood. Otherwise, PRIMA would not exist in the first place.

Thank you.

[nbelakovski]

1. Will this necessitate any change to the Fortran code?

Fortunately no.

2. For claims like the following, we need a reference to the standard.

From the 2003 standard:

Also section 13.7.88 indicates that using => null() on a pointer results in a disassociated pointer.

Link: https://wg5-fortran.org/N1601-N1650/N1601.pdf

BTW, excuse me for my ignorance about C, could you remind me why it is an arrow instead of a dot?

The arrow is used when the left side is a pointer. An example:

prima_problem_t problem;
problem.f0 = 5;
prima_problem_t *problem_ptr = &problem; // problem_ptr is of type pointer and points to problem
printf("%g\n", problem_ptr->f0); // prints 5

This is unnatural and unfriendly. If this is inevitable, then we should add another layer to enable users to do options->rhobeg=1.0. N.B.: PRIMA tries to provide the most natural and friendly interface to users and handles everything else under the hood. Otherwise, PRIMA would not exist in the first place.

OK I have an idea, the options struct could look like this (using rhobeg as an example):

typedef struct {
  double rhobeg;
  double* _rhobeg; // internal - do not use
} prima_options_t;

prima_init_options will initialize rhobeg to NAN and _rhobeg to NULL. Then the user can overwrite rhobeg and in prima_minimize we do something like if (isnan(options.rhobeg) == false) {options._rhobeg = &options.rhobeg;}. If the user hasn't overwritten rhobeg then _rhobeg will be NULL and we can set its corresponding Fortran pointer to null(). What do you think?

[zaikunzhang]

Thank you for the explanation.

This is unnatural and unfriendly. If this is inevitable, then we should add another layer to enable users to do options->rhobeg=1.0. N.B.: PRIMA tries to provide the most natural and friendly interface to users and handles everything else under the hood. Otherwise, PRIMA would not exist in the first place.

OK I have an idea, the options struct could look like this (using rhobeg as an example):

typedef struct {
  double rhobeg;
  double* _rhobeg; // internal - do not use
} prima_options_t;

prima_init_options will initialize rhobeg to NAN and _rhobeg to NULL. Then the user can overwrite rhobeg and in prima_minimize we do something like if (isnan(options.rhobeg) == false) {options._rhobeg = &options.rhobeg;}. If the user hasn't overwritten rhobeg then _rhobeg will be NULL and we can set its corresponding Fortran pointer to null(). What do you think?

Good. What about rhobeg_loc instead of _rhobeg?

Thanks.

[zaikunzhang]

I guess the following code will not compile or work as expected if RP is not C_DOUBLE.

Always keep in mind that RP is not necessarily C_DOUBLE or the default real type, and IK is not necessarily C_INT or the default integer type. Explicit conversions are always needed when we exchange data between C and Fortran (and probably between other languages). See https://github.com/libprima/prima/issues/129.

if (C_ASSOCIATED(xl)) then
    call C_F_POINTER(xl, xl_loc, shape=[n])
else
    xl_loc => null()
end if
if (C_ASSOCIATED(xu)) then
    call C_F_POINTER(xu, xu_loc, shape=[n])
else
    xu_loc => null()
end if

BTW, what does [n] mean? Will n work?

Thanks.

[nbelakovski]

I was able to make a test program work with a different value of RP by going through an intermediate step:

https://godbolt.org/z/aKdWjzYo7

Basically it works like this:

! Compulsory arguments
type(C_PTR), value :: array_of_double_ptr
integer(C_INT), value :: num_elements

! Local variables
real(C_DOUBLE), pointer :: array_of_double_ptr_loc_intermediate(:)
real(RP), allocatable :: array_of_double_ptr_loc(:)

! Convert
if (C_ASSOCIATED(array_of_double_ptr)) then
  call C_F_POINTER(array_of_double_ptr, array_of_double_ptr_loc_intermediate, shape=[num_elements])
  allocate(array_of_double_ptr_loc(num_elements))
  array_of_double_ptr_loc = real(array_of_double_ptr_loc_intermediate, kind(array_of_double_ptr_loc))
end if

There's another way to do it with pointers if necessary, but it requires a second intermediate variable. Fortunately a variable marked allocatable but not allocated will return false when present is called, just like with unassociated pointers, so we can still use this to basically use NULL to indicate than an optional variable is not provided. Should I go ahead and implement this in the PR? If so, should I go ahead and implement it for all optional variables, or just the 4 I started with?

BTW, what does [n] mean? Will n work?

From the documentation:

SHAPE | (Optional) Rank-one array of type INTEGER with INTENT(IN). It shall be present if and only if fptr is an array. The size must be equal to the rank of fptr.

Basically the brackets turn it into an array. I tried without the brackets and it complained about the variable not having the correct rank.

Good. What about rhobeg_loc instead of _rhobeg?

Usually a leading underscore is understood to mean that something is private. Plus I don't like the _loc suffix. Apparently it's short for local? That's not obvious at all, and I don't think it makes sense in this context anyway since the two are next to each other.

[zaikunzhang]

! Compulsory arguments
type(C_PTR), value :: array_of_double_ptr
integer(C_INT), value :: num_elements

! Local variables
real(C_DOUBLE), pointer :: array_of_double_ptr_loc_intermediate(:)
real(RP), allocatable :: array_of_double_ptr_loc(:)

! Convert
if (C_ASSOCIATED(array_of_double_ptr)) then
  call C_F_POINTER(array_of_double_ptr, array_of_double_ptr_loc_intermediate, shape=[num_elements])
  allocate(array_of_double_ptr_loc(num_elements))
  array_of_double_ptr_loc = real(array_of_double_ptr_loc_intermediate, kind(array_of_double_ptr_loc))
end if

This matches my imagination. Let us use safealloc instead of allocate.

Usually a leading underscore is understood to mean that something is private.

I checked this on Internet. Understood. Let us follow this convention.

Plus I don't like the _loc suffix.

I will be happy to change if there is a better option.

Thank you.

[nbelakovski]

This matches my imagination. Let us use safealloc instead of allocate.

Of course, I was only using allocate for the examples.

Should this PR make the change for all the optional variables, or just the 4 I started with? We can also make an issue to fix the rest of them in a separate PR.

[zaikunzhang]

What are the other variables?

I guess the motivation for these four variably is that they are vectors, so allocation is needed somewhere. For scalar variables, what is the motivation?

I have no objection, but I hope to understand the situation better. If the modification leads to a better (more friendly) interface or better (cleaner or more consistent) code, then of course we should do it.

It is fine to handle the other variables on a separate PR, but I hope the delay would not be long, and I hope the C API will stabilize soon.

Thank you.

[nbelakovski]

The other variables currently exposed are

• rhobeg
• rhoend
• maxfun
• iprint
• ftarget
• npt

There are of course other optional variables in the various algorithms that are not currently exposed, like honour_x0, eta, gamma, and others. For variables like maxfun and npt, we currently have logic within the C interface to initialize those, but this logic is duplicated as compared to what's in Fortran. It would be nice if C didn't duplicate any initialization code.

The motivation really isn't related to whether or not the variables are vectors or scalars. The vector quantities were initially motivated by wanting to reduce unnecessary memory allocations, but overall both types would benefit from the ability to mark them as null in order to indicate to Fortran that they are "not present."

[zaikunzhang]

For variables like maxfun and npt, we currently have logic within the C interface to initialize those, but this logic is duplicated as compared to what's in Fortran. It would be nice if C didn't duplicate any initialization code.

The motivation really isn't related to whether or not the variables are vectors or scalars. The vector quantities were initially motivated by wanting to reduce unnecessary memory allocations, but overall both types would benefit from the ability to mark them as null in order to indicate to Fortran that they are "not present."

I don't necessarily disagree. We can do it.

However, duplicating the initialization in C is not necessarily a bad thing: if a pure C implementation will be done in the future, then it will need an interface that "duplicates (implements) the initialization".

Indeed, the MATLAB interface of PRIMA and the Python interface of PDFO both implement the initialization. Indeed, they do much more, including problem cleaning, pre-processing, and refining, all of which are done very carefully. Obviously, these things are much easier to do in MATLAB/Python than in Fortran (or C). In comparison, the initialization in Fortran is relatively "primitive". I hope the Python interface for PRIMA will be similar to the MATLAB interface for PRIMA and the Python interface for PDFO.

[zaikunzhang]

I was able to make a test program work with a different value of RP by going through an intermediate step:

How to test similar things systematically? We should include such tests in our CI. See https://github.com/libprima/prima/issues/129.

[nbelakovski]

If the initialization only involved default values, then there wouldn't be very much to discuss - it would be trivial to implement in any language and not too difficult to copy the dozen or so optional variables that we have. However, I've noticed that some variables have more complicated initialization routine, like rhobeg and rhoend, for example.

	present(rhobeg) and present(rhoend)	present(rhobeg) and not present(rhoend)	not present(rhobeg) and present(rhoend)	not present(rhobeg) and not present(rhoend)
rhobeg	User-specified	User-specified	max(10*rhoend, RHOBEG_DFT)	RHOBEG_DFT
rhoend	User-specified	max(EPS, min(0.1*rhobeg, RHOEND_DFT))	User-specified	RHOEND_DFT

eta1 and eta2 have similar logic where they are initialized together. Right now in the C interface we have no way of detecting that rhoend was provided but rhobeg was not. I don't see this applying to other variables, so I think the logical course would be to implement the _rhobeg structure for these 4 variables (rhobeg, rhoend, eta1, eta2), but the others don't need this structure since they just take a default value. Since we don't currently expose eta1, eta2, I'll just implement it for rhoend and rhobeg.

[nbelakovski]

It seems nvfortran does not like my latest approach. I am investigating.

[zaikunzhang]

If the initialization only involved default values, then there wouldn't be very much to discuss - it would be trivial to implement in any language and not too difficult to copy the dozen or so optional variables that we have. However, I've noticed that some variables have more complicated initialization routine, like rhobeg and rhoend, for example.

present(rhobeg) and present(rhoend) present(rhobeg) and not present(rhoend) not present(rhobeg) and present(rhoend) not present(rhobeg) and not present(rhoend) rhobeg User-specified User-specified max(10rhoend, RHOBEG_DFT) RHOBEG_DFT rhoend User-specified max(EPS, min(0.1rhobeg, RHOEND_DFT)) User-specified RHOEND_DFT

Excuse me, @nbelakovski , it is a bit unclear to me what is the issue here.

Do you mean you will implement the same initialization in the C interface, given that we have agreed on your proposed way of indicating the absence of rhobeg (and rhoend)?

If the answer to the above question is yes, then wouldn't it contradict the following idea? According to this idea, shouldn't we just read the values of rhobeg and rhoend dumbly if they are present, and then pass whatever received to the Fortran code?

I think that ideally C should be a 'dumb' pass-through and should let Fortran handle as much as possible.

Thanks.

[nbelakovski]

Do you mean you will implement the same initialization in the C interface, given that we have agreed on your proposed way of indicating the absence of rhobeg (and rhoend)?

No, I was pointing out that the C interface as currently designed (i.e in the main branch, not in this PR) is broken in the sense that it can only pass values for rhobeg and rhoend and cannot indicate if they are absent. So if the user supplies rhobeg and does not supply rhoend, instead of setting rhoend to 0.1*rhobeg, it will just use the default value for rhoend as currently specified in the C interface.

This seems like a problem to me, and given my statement that I think C should be a 'dumb' pass-through, I think the way to solve this is to use that approach with rhobeg and _rhobeg to indicate to Fortran whether the user has provided a value for that variable or not. Then Fortran can do the proper initialization. Does that make sense?

I've investigated the nvfortran issue. The issue seems to come from the -Mchkptr flag passed to nvfortran. It complains when bobyqa is called from bobyqa_c that rhobeg_loc is null. This seems like a bug with nvfortran, since we're not dereferencing the pointer at this stage (and in fact in Fortran it's legal to pass unitialized or unallocated pointers to intrinsic procedures like present). I don't see a way around this in code. Do you think it would be appropriate to either remove the -Mchkptr flag from cmake.yml, or to modify the CMake files so that the flag is used when compiling Fortran examples but not the C interface?

[zaikunzhang]

if the user supplies rhobeg and does not supply rhoend, instead of setting rhoend to 0.1*rhobeg, it will just use the default value for rhoend as currently specified in the C interface.

You are right. This setting is wrong.

This seems like a problem to me, and given my statement that I think C should be a 'dumb' pass-through, I think the way to solve this is to use that approach with rhobeg and _rhobeg to indicate to Fortran whether the user has provided a value for that variable or not. Then Fortran can do the proper initialization. Does that make sense?

Yes. Agreed.

[zaikunzhang]

Thank you @nbelakovski for the investigation.

This seems like a bug with nvfortran, since we're not dereferencing the pointer at this stage (and in fact in Fortran it's legal to pass unitialized or unallocated pointers to intrinsic procedures like present).

I see. It will not surprise me if we have caught another compiler bug. It would be better if we can submit a minimal example to nvfortran and get a confirmation that it is their bug.

Do you think it would be appropriate to either remove the -Mchkptr flag from cmake.yml, or to modify the CMake files so that the flag is used when compiling Fortran examples but not the C interface?

The second approach sounds better unless it is much more complicated to implement.

[nbelakovski]

It looks like someone has run into this before: https://forums.developer.nvidia.com/t/bug-in-nvfortran-with-mchkptr-for-unallocated-optional-arguments/223220

They claim it's correct behavior, but they understand that this is a valid use case and offered to make a second flag to allow this sort of behavior while -Mchkptr is used. It doesn't seem like there was any follow up. I'll make a topic to follow up.

In the meantime, I dug around, and it seems like it might technically be possible to remove a flag from within CMake for a specific target, but there's no standard-conforming way to do it, and I think it would be best to just remove the chkptr flag from CMake. If we need it, we can make a separate github action that only builds and tests the Fortran code and not the C interface.

Edit: nvidia forum topic to revive the discussion around Mchkptr suboption https://forums.developer.nvidia.com/t/mchkptr-doesnt-allow-passing-null-optional-values/278293

[nbelakovski]

OK this PR is good to go @zaikunzhang

[zaikunzhang]

Thank you @nbelakovski . I have made some comments. Can you see them (I am not sure)?

Do not do a force push before checking all the comments, or they will disappear.

[nbelakovski]

No, I didn't see them, I'm sorry I didn't know force-pushing removed comments. Could you please re-comment?

[zaikunzhang]

What about now? @nbelakovski

[nbelakovski]

@zaikunzhang I've pushed a new version

• _intermediate was changed to intrmdiate to try to match the style of shortened variable names
• rhobeg and rhoend and f0 were modified to have nan as their default value in C, and are marked as allocatable in Fortran so that if the C variable is nan, then rhobeg_loc and rhoend_loc will be unallocated and will therefore be interpreted as not present. I thought this was a reasonable approach since your comment in the 0.7.2 release indicated that is_nan resolved the issues with gcc 13 -Ofast.

[zaikunzhang]

_intermediate was changed to intrmdiate to try to match the style of shortened variable names

Sorry, I don't understand this change. I guess most people will interpret intrmdiate as a typo. It is worse than before.

See https://www.google.com.hk/search?q=intermedite+abbreviation&ie=UTF-8&oe=UTF-8&hl=fr-hk&client=safari#sbfbu=1&pi=intermediate%20abbreviation

[zaikunzhang]

It seems that some comments are deleted. There were some points I have not responded to yet.

[zaikunzhang]

It seems that some comments have been deleted. There were some points I have not responded to yet.

Let me respond according to my memory:

What if rhobeg = NaN, but is_nan does not work, so that rhobeg is then set to a huge value?

Answer: It will be a disaster. The algorithm will still work in theory, but it will take a loooong time to converge. In real computation, this essentially means that the algorithm will not do anything correct.

[zaikunzhang]

is_nan resolved the issues with gcc 13 -Ofast.

There are still cases where is_nan does not work. For example, https://github.com/equipez/infnan/issues/29.

[nbelakovski]

It seems that some comments are deleted. There were some points I have not responded to yet.

They are not deleted, I've simply marked some of them as resolved. You can open them up again by clicking "Show resolved"

[nbelakovski]

I'm a little confused as to why it is difficult to check for NaN? From wiki, the IEEE 754 standard indicates NaN by setting the exponent to all 1's and the significand to a nonzero value. Isn't it possible to do a bitwise and with a number with all 1's in both the sigificand and the exponent and check that the resulting value is nonzero?

[zaikunzhang]

I'm a little confused as to why it is difficult to check for NaN? From wiki, the IEEE 754 standard indicates NaN by setting the exponent to all 1's and the significand to a nonzero value. Isn't it possible to do a bitwise and with a number with all 1's in both the sigificand and the exponent and check that the resulting value is nonzero?

At the beginning, I was surprised as well. Now nothing like this can surprise me anymore.

Fortran standard says nothing about IEEE standards. Compilers do not necessarily respect IEEE standards, and they actually don't --- this is what I have learned in the past years.

Recall that we have three real kinds to handle. This adds more complexity to the problem.

Isn't it possible to do a bitwise and with a number with all 1's in both the sigificand and the exponent and check that the resulting value is nonzero?

There is no way to record such a number if the real is quadruple.

[nbelakovski]

I'm a little confused as to why it is difficult to check for NaN? From wiki, the IEEE 754 standard indicates NaN by setting the exponent to all 1's and the significand to a nonzero value. Isn't it possible to do a bitwise and with a number with all 1's in both the sigificand and the exponent and check that the resulting value is nonzero?

At the beginning, I was surprised as well. Now nothing like this can surprise me anymore.

Fortran standard says nothing about IEEE standards. Compilers do not necessarily respect IEEE standards, and they actually don't --- this is what I have learned in the past years.

Recalling that we have three real kinds to handle. This adds more complexity to the problem.

In this case the NAN is being defined in C and then passed to Fortran, does that change things?

[nbelakovski]

_intermediate was changed to intrmdiate to try to match the style of shortened variable names

Sorry, I don't understand this change. I guess most people will interpret intrmdiate as a typo. It is worse than before.

See https://www.google.com.hk/search?q=intermedite+abbreviation&ie=UTF-8&oe=UTF-8&hl=fr-hk&client=safari#sbfbu=1&pi=intermediate%20abbreviation

None of those abbreviations make any sense. I don't find the style argument very compelling either. Why make something unintelligible just because the things nearby are unintelligible? If that was the case you could have stayed with the original F77 code :joy:

[nbelakovski]

I've explored the nan issue and I see that it's quite messy and that it's an issue in C as well as Fortran.

In light of this I've made the following changes:

• Since rhobeg and rhoend need to be positive, I've set their initialize value at -1 and in algorithm_c.f90 I check if the value of the option is >0 and then I set it, otherwise it's left unallocated and will be interpreted as not present.
• I left the NaN functionality for f0 since it can take on any value. My rationale is that the NaN approach leads to friendlier code for the user (they can set problem.f0 = ... instead of double f0 = ...; problem.f0 = &f0;) and in exploring COBYLA I got the impression that the way it tracks multiple values in a simplex might make it robust in the case of f0 accidentally being a huge value? Of course I've only tested it on a trivial problem so I'm looking to you for guidance here. If you believe it's not robust to f0 accidentally being huge then we can go back to making f0 a pointer (this probably won't affect the interface much since I think most users will initialize it by passing the address of problem.f0 to the objective function, as opposed to using a separate value to initialize f0).
• I changed the variable name to interm because I don't want to get stuck on this.

[zaikunzhang]

Thank you @nbelakovski .

if (.not. is_nan(f0)) then
    f0_loc = real(f0, kind(f0_loc))
end if
...
if (rhobeg > 0) then
    rhobeg_loc = real(rhobeg, kind(rhobeg_loc))
end if
if (rhoend > 0) then
    rhoend_loc = real(rhoend, kind(rhoend_loc))
end if

Let's not do this.

1. The only valid way of indicating missing data is to use something that means "missing" or "unset", e.g., NaN (for number), NULL (for pointer), none (Python), or simply non-existence of the variable (doable in Python and Matlab, not in C or Fortran). It is logically WRONG to indicate missing data by invalid data. It may be fine for the moment, but it will lead to maintenance difficulty in the future after some code evolution. I know it.

2. We should not assign unnatural meaning to data. By unnatural, I mean anything artificial and nonessential. rhobeg = -1 should mean nothing but setting rhobeg to minus one, which should only be interpreted as nothing but an invalid input.

3. The Fortran code treats rhobeg = -1 as an invalid input, warns about it, and corrects it. See https://github.com/libprima/prima/blob/f92eb14d94f7d46eaf18a661d07d39332382efaa/fortran/common/preproc.f90#L297-L306 The C interface should either do the same or do nothing. Otherwise, the code will quickly become non-understandable.

4. Moreover, In this piece of code, we are using two different ways to indicate missing data (NaN and an invalid value). This is not good.

5. What will happen if f0 is mistakenly set to a huge value? Answer: Disaster. The method is based on interpolation of the function values. Any small mistake in the function value will hurt the performance. A huge mistake will kill the algorithm. (Of course, in real problems, the function values often contain noise --- this is something inevitable and we have to tolerate, but mistakenly setting f0 to a huge value is different).

6. If we had an isnan that always works, then I would prefer setting xyz = NaN to indicate that xyz is not provided. Otherwise, I think your previous proposal is good (with some comments):

typedef struct {
  double rhobeg;
  double* _rhobeg; // internal - do not use
} prima_options_t;

I changed the variable name to interm because I don't want to get stuck on this.

Thank you.

[nbelakovski]

It is logically WRONG to indicate missing data by invalid data. It may be fine for the moment, but it will lead to maintenance difficulty in the future after some code evolution. I know it.

I'm not convinced of this. In cases where data is supposed to be positive, using negative to indicate not available seems like a reasonably solution if other options are not available. PRIMA is already doing this with maxfun and npt. However this is not the first time I've had this conversation and I've come to realize that some people hate this approach and some people think it's fine. Clearly you're in the former camp so I won't try to convince you otherwise.

If we had an isnan that always works, then I would prefer setting xyz = NaN to indicate that xyz is not provided. Otherwise, I think your previous proposal is good (with some comments):

I don't see any comments, it looks like you just reposted my previous code? More importantly there remains an issue, the approach with rhobeg and _rhobeg still relies on NaN. It sets rhobeg to NaN and _rhobeg to NULL, and only sets _rhobeg to a valid address if rhobeg is not NaN. The only way to guarantee using NULL is to ask the user to manage the memory (i.e. by doing double rhobeg = 1; options.rhobeg = &rhobeg;). It's not clear to me where your priority lies in terms of user friendliness and correct operation across all compiler options, and so I don't know which option to implement. I did find a way to detect -ffast-math in C code (see here), but I don't know which C compilers and versions you'd like to support, so again I don't know if I should use isnan and check for fast-math, or if I should assume that there is some compiler+option out there that will break NaN and therefore should avoid its use entirely.

So in summary I need some guidance since I don't know how solve all of the problems you have posed, I can only solve some of them and I'm not sure which ones to solve.

[zaikunzhang]

It is logically WRONG to indicate missing data by invalid data. It may be fine for the moment, but it will lead to maintenance difficulty in the future after some code evolution. I know it.

I'm not convinced of this. In cases where data is supposed to be positive, using negative to indicate not available seems like a reasonably solution if other options are not available.

Sorry, this approach is wrong. People are doing something wrong just means they are making mistakes. It doesn't justify the wrong approach. But as you said, I will not try to convince you.

PRIMA is already doing this with maxfun and npt.

No. I suppose you mean the following code. If you read it again, you will say that it treats a negative value as an invalid input, warns about it (note the warning), and corrects it. This is different from treating a negative value as a void/non-existent input.

if (maxfun <= max(0, min_maxfun)) then
    if (maxfun > 0) then
        maxfun = int(min_maxfun, kind(maxfun))
    else  ! We assume that nonpositive values of MAXFUN are produced by overflow.
        if (MAXFUN_DIM_DFT >= huge(maxfun) / n) then  ! Avoid overflow when N is large.
            maxfun = huge(maxfun)
        else
            maxfun = MAXFUN_DIM_DFT * n
        end if
    end if
    call warning(solver, 'Invalid MAXFUN; it should be at least '//min_maxfun_str//'; it is set to '//num2str(maxfun))
end if

! Validate NPT
if ((lower(solver) == 'newuoa' .or. lower(solver) == 'bobyqa' .or. lower(solver) == 'lincoa') &
    & .and. present(npt)) then
    if (npt < n + 2 .or. npt >= maxfun .or. 2 * npt > int(n + 2) * int(n + 1)) then  ! INT(*) avoids overflow when IK is 16-bit.
        npt = int(min(maxfun - 1, 2 * n + 1), kind(npt))
        call warning(solver, 'Invalid NPT; it should be an integer in the interval [N+2, (N+1)(N+2)/2]'// &
            & ' and less than MAXFUN; it is set to '//num2str(npt))
    end if
end if

However this is not the first time I've had this conversation and I've come to realize that some people hate this approach and some people think it's fine. Clearly you're in the former camp so I won't try to convince you otherwise.

Agree. I will do the same.

If we had an isnan that always works, then I would prefer setting xyz = NaN to indicate that xyz is not provided. Otherwise, I think your previous proposal is good (with some comments): ' I don't see any comments, it looks like you just reposted my previous code?

Yes, I copy-pasted your old code.

I think your previous proposal is good (with some comments):

My English was not clear. I meant this: your previous approach is good as long as it is accompanied with sufficient comments.

More importantly there remains an issue, the approach with rhobeg and _rhobeg still relies on NaN. It sets rhobeg to NaN and _rhobeg to NULL, and only sets _rhobeg to a valid address if rhobeg is not NaN. The only way to guarantee using NULL is to ask the user to manage the memory (i.e. by doing double rhobeg = 1; options.rhobeg = &rhobeg;).

I am confused. You mentioned that we can pass a null pointer to Fortran, which will be interpreted as absence of the argument. Where do we need the NaN?

It's not clear to me where your priority lies in terms of user friendliness and correct operation across all compiler options, and so I don't know which option to implement.

Maximise the portability subject to user friendliness and code quality&clarity. User friendliness and code quality&clarity are the priority and they are not contradicting each other.

I have a very strong preference for user friendliness and code quality&clarity because I have lost too much time due to some codebase that sacrificed them for something that are nonessential from a modern point of view. If user friendliness and code quality&clarity were not our priority, PRIMA would not exist in the first place.

If some compiler options cannot be supported under these constraints due to compiler bugs/limitations, we can ignore them as long as we are sure it is not a bug/limitation of our code.

So in summary I need some guidance since I don't know how solve all of the problems you have posed, I can only solve some of them and I'm not sure which ones to solve.

If my understanding is correct, your previous approach would solve all problems.

Thank you.

[nbelakovski]

No. I suppose you mean the following code.

No, I meant the C interface code.

I am confused. You mentioned that we can pass a null pointer to Fortran, which will be interpreted as absence of the argument. Where do we need the NaN?

The NaN is not in Fortran, it is in C. C will initialize rhobeg to NaN and _rhobeg to NULL. _rhobeg is what gets passed to Fortran, so if it is null it will be interpreted as absent. But in this situation how does the user provide rhobeg? They setrhobeg to some desired value, and when they called prima_minimize, the code will check if rhobeg is NaN (the default value), and if it is not NaN, then it will set _rhobeg = &(options.rhobeg);. So in C we need to check if rhobeg is NaN. Fortran will only see a pointer. But if the code is compiled with gcc with -ffast-math this check will fail and _rhobeg will end up pointing to a large value. Does that make sense? clang appears to preserve isnan even with fast-math and I was able to find a way to detect when the code is compiled with gcc and fast-math, but like I said earlier I don't know if there are other C compilers I need to test, or what to do if I detect gcc and fast-math (error? warning? custom isnan implementation?).

[zaikunzhang]

No. I suppose you mean the following code.

No, I meant the C interface code.

Which part? Then it is a mistake and should/must be corrected.

[nbelakovski]

No. I suppose you mean the following code.

No, I meant the C interface code.

Which part? Then it is a mistake and should/must be corrected.

https://github.com/libprima/prima/blob/main/c/prima.c#L17 and https://github.com/libprima/prima/blob/main/c/prima.c#L102, but what about my comment regarding NaN and detecting fast-math?

[zaikunzhang]

I was able to find a way to detect when the code is compiled with gcc and fast-math.

Let us not do this kind of detection.

On the Fortran side, the code is tested with 4 optimization levels. For each level, it is tested that the code compiles without any warning and works correctly.

-g: corresponding to the lowest optimization level of the compiler, nomally compiler -g -O0 -O1: corresponding to compiler -O1 -O2: corresponding to compiler -O2 -O3: corresponding to compiler -O3 -fast: corresponding to the highest optimization level of the compiler, e.g., gfortran -Ofast, ifx -fast, etc

Which C compilers to test? I would test as many as possible. This includes at least gcc, icc, and clang. Are there others? I do not know.

[zaikunzhang]

OK. I comprise. Let us use NaN to indicate an absent real scalar , Null to indicate an absent vector (what about something like [] or {}?). What is your suggestion for absent integers? @nbelakovski

[nbelakovski]

Let us not do this kind of detection.

May I ask why not? I'd be reluctant to make code modifications with this kind of detection, as it quietly defeats the purpose of these flags, but I think that a warning or an error might be appropriate for this sort of circumstance (particularly since many people complain about gcc -Ofast breaking isnan, in fact so many people complained about it that clang decided that -Ofast would not break isnan). I'm not really trying to push for warning or error here, more just interested in your thoughts.

What is your suggestion for absent integers?

Since integers can't be NaN, I think the only option is NULL. In that case there would be no way to avoid forcing C users to manage the memory for things like npt, should they choose to specify it.

what about something like [] or {}?

double* vect = {}; will lead vect to be null, but in our case we memset the whole problems and options structs to 0 so that sets any pointer to null.

Edit: I've now pushed changes to reflect these discussions. I wasn't entirely sure if you preferred the is_nan check to be in Fortran or in C so I put it in Fortran so that we don't have to use the _rhobeg method in C and complicate the code, but as always I'm open to reinterpretation if you wish. With that, do we have any open questions left on this PR?

[zaikunzhang]

Let us not do this kind of detection.

May I ask why not? I'd be reluctant to make code modifications with this kind of detection, as it quietly defeats the purpose of these flags, but I think that a warning or an error might be appropriate for this sort of circumstance (particularly since many people complain about gcc -Ofast breaking isnan, in fact so many people complained about it that clang decided that -Ofast would not break isnan). I'm not really trying to push for warning or error here, more just interested in your thoughts.

OK. A warning in the building system (now CMake, later Meson) would be appropriate. But, as you said, we will not change the C/Fortran code according to the compiler flags. The warning comes from the building system, not the solvers.

What about "Warning: compiling the library with gcc -Ofast; isnan may not work and the solvers may behave incorrectly."

This is needed only if we check NaN in C, right (see below).

What is your suggestion for absent integers?

Since integers can't be NaN, I think the only option is NULL. In that case there would be no way to avoid forcing C users to manage the memory for things like npt, should they choose to specify it.

Let us use -1 for npt, maxfun, maxhist, and maxfilt :smile: , the last two being not exposed for the moment.

Except for iprint, all integer arguments of PRIMA are expected to be positive integers. Let us reserve -1 (and only -1) to indicate their absence; other nonpositive values will be passed to Fortran as is, and will be warned and corrected.

As for iprint, it takes 0, +-1, +-2, and +-3; any value will be passed to Fortran as is (it will be warned and corrected if invalid).

what about something like [] or {}?

double* vect = {}; will lead vect to be null, but in our case we memset the whole problems and options structs to 0 so that sets any pointer to null.

OK. You will choose the implementation that is the best according to your expertise.

Edit: I've now pushed changes to reflect these discussions. I wasn't entirely sure if you preferred the is_nan check to be in Fortran or in C so I put it in Fortran so that we don't have to use the _rhobeg method in C and complicate the code, but as always I'm open to reinterpretation if you wish. With that, do we have any open questions left on this PR?

I have not checked the code. I suppose you meant the bind-C Fortran code (solve_c) when you said "Fortran", right? Then yes, I think it is the correct place. We should use the is_nan function from infnan_mod. In this way, we do not need to worry about the isnan in C. Consequently, we do not need a warning about gcc -Ofast.

What about a warning about gfortran -Ofast for the Fortran code? I do not think it is necessary for the moment. The Fortran is_nan function is robust enough. When it does not work, I tend to regard it as a bug of the compiler, for which we do not assume responsibility.

Thank you.

[zaikunzhang]

Let us use -1 for npt, maxfun, maxhist, and maxfilt 😄 , the last two being not exposed for the moment.

Except for iprint, all integer arguments of PRIMA are expected to be positive integers. Let us reserve -1 (and only -1) to indicate their absence; other nonpositive values will be passed to Fortran as is, and will be warned and corrected.

As for iprint, it takes 0, +-1, +-2, and +-3; any value will be passed to Fortran as is (it will be warned and corrected if invalid).

@nbelakovski

I have a better idea. For all these optional integer arguments, we use 0 to indicate that the argument is absent & should take the default value.

For optional real arguments, we use NaN.

For arrays, we use NULL.

I will change the Fortran backend to implement the first two points. The third one does not need changes to the Fortran backend. The Fortran backend will be silent when receiving, e.g., npt = 0 or rhobeg = NaN.

After these changes, the C interface can simply pass whatever is received to Fortran.

The MATLAB/Python/Julia/R interface, however, should not use 0 or NaN to indicate the absence of an option. They should (eventually) have a complete preprocessing procedure that interprets and preprocesses the user's input (see preprima.m, especially the preoptions function). 0 and NaN will be interpreted as invalid values, warned, and corrected. The preprocessing procedure will need some time to code. We can first start with something simple but works. preprima.m is the reference implementation for such a procedure.

In other words, we use 0 and NaN to indicate absence only in the C interface, due to the limitation of the language. Other languages have better ways to do the same thing and do not need this compromise.

[zaikunzhang]

I will change the Fortran backend to implement the first two points. The third one does not need changes to the Fortran backend. The Fortran backend will be silent when receiving, e.g., npt = 0 or rhobeg = NaN.

Hi @nbelakovski , I have done this: https://github.com/libprima/prima/commit/8eb094ad3d3f6497b8f6821d96bccfd75b788b07 , https://github.com/libprima/prima/commit/ef6e31788ec6b50e01f043252aa7ede87f86664c

[nbelakovski]

@zaikunzhang I've made the relevant changes as well. I've also added some notes in c/README.md regarding our discussion of default values and linked to this PR.

Is there anything else?

[zaikunzhang]

I do not think it is necessary to check whether rhobeg or rhoend is NaN, given

https://github.com/libprima/prima/blob/f55179c0a4ef702eccd3e6f90637168e346ca044/fortran/common/preproc.f90#L312-L331

For any real scalar, the Fortran code will interpret NaN as an absent argument. Hence the C interface can simply pass whatever is received to the Fortran backend.

Recall #135 (comment) ("After these changes, the C interface can simply pass whatever is received to Fortran.")

What do you think?

Thank you.

On my side, this is the only thing left before merging.

[nbelakovski]

On my side, this is the only thing left before merging.

For this I'd apply the same arguments I applied towards making sure we initialize npt and maxfun within the C code. The reason to use NaN to indicate that they're absent is so that the Fortran code can take the default path, instead of passing NaN through the whole chain of prima_minimize -> algorithm_c -> algorithm -> preproc before it gets acted upon.

[zaikunzhang]

On my side, this is the only thing left before merging.

For this I'd apply the same arguments I applied towards making sure we initialize npt and maxfun within the C code. The reason to use NaN to indicate that they're absent is so that the Fortran code can take the default path, instead of passing NaN through the whole chain of prima_minimize -> algorithm_c -> algorithm -> preproc before it gets acted upon.

But this substantially complicates the code. In addition, the code about f0 is wrong - check the code about f0 and nlconstr0 in cobyla.f90 (I am on cellphone so I cannot make a link easily). Thanks.