I’ve been using stdlib in new projects and replaced older linear algebra procedures (my own old poorly optimised ones or standalone lapack) with the stdlib_linalg interfaces as much as possible. I like to keep the core calculations in pure procedures. However, I rely quite a bit on eigendecomposition and stdlib’s eigh, for example, is impure.
My questions:
For those who did all the impressive work on stdlib_linalg: What determined whether procedures were made pure (e.g., cholesky) or left impure (as eigh and equivalents in the old lapack)? (I’m also trying to get a sense of expected challenges if I decide to dive into that part of the code).
More general: Does anyone have suggestions for a workaround (e.g., similarly well written, MIT licenced pure procedures)?
When converting LAPACK to Modern Fortran, the pure and elemental attributions were made from the bottom up, i.e., starting from the simplest “auxiliary” functions and then propagating the pure attribute to all procedures using them.
Unfortunately the culprit that breaks total purity is the original LAPACK interface containing several functions that have intent(out) or intent(inout) arguments.
The top-level LAPACK interface cannot be changed, but a workaround would be to make said functions just wrappers over pure subroutines, and then changing all places they’re called in the code with calls to the subroutines instead of the functions.
Think something like
function original_interface(a,b,c) result(d)
real, intent(inout) :: a(*),b(*),c(*) ! not pure
real :: d
call original_interface_sub(a,b,c,d)
end function
pure subroutine original_interface_sub(a,b,c,d)
real, intent(inout) :: a(*),b(*),c(*) ! can now be pure
real, intent(out) :: d
end subroutine original_interface_sub
I believe the reason for this constraint is that functions are “special”, i.e., the compiler may optimize the call order of say, f(x,y) + f(x,y+1): if x is modified, the order of the calls matter.
Such an expression is disallowed by 10.1.4 Evaluation of operations, para 2 (J3/25-007r1), which includes:
• the evaluation of a function reference shall neither affect nor be affected by the evaluation of any other entity within the statement;
• if a function reference causes definition or undefinition of an actual argument of the function, that argument or any associated entities shall not appear elsewhere in the same statement.
If that is forbidden anyways, why then enforce intent(in) arguments in functions? Could not functions have constant references to non-constant data via intent(inout) like subroutines can? the no-aliasing rule always applies anyways in Fortran
My question is: if the only case where intent(inout) would bother a pure function (see the f(x,y) + f(x,y+1) example with intent(inout) :: x`) is also already forbidden by another rule, why enforce it further?
The reason I’m asking is, it’s natural and native to think about functions as math functions in Fortran, but let’s not forget that other languages often use function results as error handlers, etc. so having a pure function doing some manipulations and returning intent(inout) arguments would have a lot of applications:
The prohibition I pointed to predates PURE and INTENT and is already present in Fortran 66, but judging from some codes I’ve seen and heard about, was seldomly (never?) enforced by compilers. A Fauxrtran “program” that violates it but is claimed by its author to “work” is a maintenance nightmare. The newer syntax was intended to enable compile-time checks to avoid getting into this situation in the first place.
The PURE attribute is supposed to do other things, such as allow constant expression values to be extracted from within loops or the simultaneous evaluation of functions in multiple threads or execution streams.
I used to follow a similar approach as you and tried to make all routines central to my calculations pure, in large part so I could use Do Concurrent loops.
However I ended up with the conclusion that pure procedures in their current form are just too limiting to be a hard requirement, and as a result I dumped Do Concurrent in all my codes and replaced them with OpenMP calls and haven’t looked back. I will label a pure if I can but it’s not a requirement.
In addition to current struggles you are working though, there are also issues with any procedure that allocates a variable not being allowed to be pure that is a struggle as well.
I haven’t found it too limiting or a problem so far (except for what’s mentioned in this thread), and I find that separating core calculations from other functionality gives it a more understandable structure and more flexibility even without any consideration given to parallelisation. However, I appreciate that the degree of (in)convenience of such an approach depends on the nature of the project.
What do you mean by the last sentence? I’m not sure I understand. A problem created by a called pure procedure allocating and returning a variable? How so?
This is discussed in several threads on the Intel forum. Looking back the restrictions on allocate in pure procedures are specific to polymorphic objects. The compiler can’t know at compile time of the object has a finalized and if it is pure. It was published as a ‘corrigendum’ to F08 I believe.
Wow, nice! Great to hear that such a systematic approach was taken to modernise the procedures in that way.