Introduce initialization and free of abstract vector types

[gloopydoop]

This PR introduces an overwritable init and free procedure for the abstract vector type.
The intention is for abstract vectors which may require sophisticated initialization and destruction, and is primarily focused towards supporting GPUs and is specifically related to the codes neko and neko-top. The following is a list of the major changes:

1. Many complex objects in neko require initialization and freeing subroutines (I think this is quite common in general). This is particularly important for our code when we run on GPUs. For us, a field is a rather complex object. It contains standard arrays when running on CPUs but when one runs on GPUs it contains pointers to the memory on device, functionality to copy between host and device etc. So the first major edit was to introduce an init and a free to the abstract_vector. By default this does nothing, and so the user doesn't have to prescribe two extra subroutines. However, these subroutines can be overridden if the user desires.

2. An unfortunate consequence of allocating with allocate(X, source=Y) is that the device pointers are copied across. So this would mean X and Y technically share the same data on the GPU. So most of these were replaced with allocate(X, mold=Y); call X%init(); call copy(X, Y).

3. A handful of X = Y were replaced by copy(X, Y).

4. This was maybe the most intrusive... There were a number of subroutines, ie, copy, which had intent(out) for X. In principal, this should give us a fresh X, so it would require an X%free(); X%init() on entry. However, since this was used so extensively I felt it would be very costly to constantly break down and rebuild device pointers for every copy so I changed to intent(inout).

For a test usage case please see ExtremeFLOW/neko-top#195.

[codecov]

Codecov Report

❌ Patch coverage is 91.52299% with 59 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/AbstractTypes/AbstractLinops.f90	55.55%	16 Missing ⚠️
src/AbstractTypes/AbstractVectors.f90	94.66%	16 Missing ⚠️
src/IterativeSolvers/CG/CG.f90	73.21%	15 Missing ⚠️
src/Krylov/gram_schmidt.f90	50.00%	12 Missing ⚠️

📢 Thoughts on this report? Let us know!

[Simkern]

Hej @gloopydoop!

I went through your PR and the suggestions you have made and I have created new branch with my proposition that essentially follows your proposition, but streamlines some of the implementation details.

The current implementation passes all our CPU tests. I do not have access to a GPU so I do not know if something breaks if we actually hook it up to a GPU. Maybe @gloopydoop can try this with neko?

Check out the pr-neko_fix_sk branch if you'd like to compare.

I'm not claiming that my implementation is better, especially since it is more intrusive that @gloopydoop 's, but I put it out there as a basis for discussion on how to proceed in this regard.

@loiseaujc : You should take a look and add your 5 cents so that we can move forward.

The changes and logic of the alterations for my version are the following:

pr-neko_fix_sk

This PR reworks the abstract vector lifecycle so generic Krylov algorithms run on CPU and GPU backends. Moves all allocation off source= (which shallow-copies device handles) onto an explicit mold= + init_like shaping step, with matching free cleanup.

New API

• init_like(self, mold): new deferred type-bound procedure on abstract_vector_{r,c}{sp,dp}. Sole shape primitive: allocate if unallocated / no-op if conformant / reallocate if shape differs. Must be idempotent.
• free(self): new overrideable type-bound procedure, default no-op (clears is_initialized). Terminal + idempotent.
• init_like_basis(X, mold), free_basis(X) — elemental basis wrappers (now public).
• abstract_vector base gains inert is_initialized / owns_data flags (scaffolding for GPU; unused + untested for CPU).

Internal changes (no user action)

• All solvers (gmres, fgmres, cg, eighs, eigs, svds, kexpm, Newton) switched source= -> mold= + init_like, with explicit free/free_basis cleanup.
• Newton: sys%jacobian%X = X -> guarded allocate(mold) once + copy in-loop.

Performance improvements (not directly related to the new API)

• kexpm_mat/kexpm_vec: projection hoisted out of the Arnoldi loop (computed once at convergence, not every step).
• permcols_basis: in-place copy(Q, Q(perm)) replaced with a temp-buffer permutation (aliasing/GPU-safe).

Bug fix

• dense_axpby with beta == 0 now overwrites (self = 0 then axpy) instead of skipping the scale. The old if (beta /= 0) scal left stale contents in self, corrupting copy into reused storage (

This bug was invisible with the previous implementation but silently broke the GMRES test with the new API.

BREAKING CHANGES

1. All user-defined vector types must implement the deferred init_like.
2. copy(out, from) now requires out to be pre-allocated with intent(inout)
3. Operator/solver output args changed intent(out) -> intent(inout) across the abstract interfaces (matvec/rmatvec/apply, system response, eigensolver/SVD X/U/V, kexpm C, krylov_exptA vec_out).

Comments

• I chose to implement a deferred type-bound procedure for init_like instead of the alternative to use the generic subroutine init_like with the select_type construct. This seems very hard to maintain and less clear than a TBP associated with each concrete type allowing for the initalization knowing all the specifics of the concrete type in question.

• @gloopydoop 's suggestion to solve the initialization issue for the matvec operation by placing the init_like call in the internal apply_matvec routine is minimally intrusive. However this is dangerous since it hides this change from the user who would normally never look at these utility wrappers that exist only to allow for execution timing. Since the allocate -> init_like chain is explicit for the vector types, it should also be explicit for the matvecs.

• The question of the init_like call site also raises a question regarding performance. If we maintain the intent(out) logic that we had initially chosen for clarity and to avoid the old fortran issues of mistaking inputs and outputs, we are always required to deallocate/reallocate at each call and cannot take advantage of a global scratch space with a one-time allocation and final free. I think this choice was defensible in the pure CPU setup since all of the allocation and memory management was essentially outsourced, it becomes harder to motivate if we include the possibility of fine grained memory control that we assume the user has anyway.

• Side note: We are now generally inconsistent with the intents. Some functions are intent(out), most prominently the linear_combination but most other subroutines are now intent(inout). This makes the code in general less clean, in particular because we have no in-place linear_combination, i.e. no wrapper that can act on an intent(inout) target. This means that we are duplicating the exact code from linear_combination in these places which is much less readable because a normal user will assume that the code is either a duplicate where we forgot to use linear_combination or that it's genuinely doing something different. I would advocate for one of two options:

  1. add a utility function for this case, e.g. linear_combination_in_place or similar to be used in these situations.
  2. fully switch to intent(inout) globally, fully transferring the memory management to the user as is often the case in Fortran.

I'm looking forward to your suggestions regarding the implementation!