linalg eye: allow generalized return type and kind

[perazz]

It was recently suggested (#901) that it may be confusing that eye always returns an integer(int8) result.
This PR proposes to extend the current identity matrix implementation via an optional mold keyword.
If present, mold is used to give a type and kind to the return matrix.

Proposed interface

• A = eye(dim1 [, dim2] [, mold])

Key facts

• Default behavior is unchanged (implementation is backward compatible)
• Fortran Discourse thread

Prior art

mold already exists, and with the same meaning, both in the Fortran standard (allocate(a, mold=b)), and here in stdlib:

• codata constants
• string to number conversion

cc: @jvdp1 @jalvesz @fortran-lang/stdlib

[jvdp1]

Thank you @perazz . This is a nice extension. I just wonder if we should keep the default output as integer(int8) or as integer?

[perazz]

@jvdp1 let me know if you have further comments since I switched the default return to real(real64). Otherwise, I believe this can be merged soon.

[jalvesz]

Thanks @perazz for this improvement! I think it is indeed almost ready. One suggestion, feel free to neglect it:

Quite often it is required to initialize arrays with zero or one, why not add these parameters to one of the constants modules? taking for instance the one added in the stdlib_sparse_constants:

#:for k1, t1, s1 in (R_KINDS_TYPES)
    ${t1}$, parameter :: zero_${s1}$ = 0._${k1}$
    ${t1}$, parameter :: one_${s1}$ = 1._${k1}$
#:endfor
#:for k1, t1, s1 in (C_KINDS_TYPES)
    ${t1}$, parameter :: zero_${s1}$ = (0._${k1}$,0._${k1}$)
    ${t1}$, parameter :: one_${s1}$ = (1._${k1}$,1._${k1}$)
#:endfor

Make them accesible elsewhere, in that case, I would propose the following small modification to the eye template:

...
allocate(result(dim1, dim2_), source = zero_${s1}$)
        
do i = 1, min(dim1, dim2_)
      result(i, i) = one_${s1}$
end do
...

Note: I approve it even if these suggestions are not used as I think it is ready as is ;)

[perazz]

Thank you @jalvesz @jvdp1, with two approvals, I will merge this one.
Regarding your suggestion: I think it's a good idea and we should think of a good way to do it (perhaps, inside the stdlib_linalg_constants module?)

[jvdp1]

Quite often it is required to initialize arrays with zero or one, why not add these parameters to one of the constants modules?

What would be the advantage of this approach for eye? Would it not complicate the code (with all the different possibilities among int, real, and complx)?
In any case, I think it should be part of another PR.

[jalvesz]

@jvdp1 I think that if consistently used, it would avoid the small kind conversion errors under assignment.

Agree, let's postpone that for another PR.

[jvdp1]

@jvdp1 I think that if consistently used, it would avoid the small kind conversion errors under assignment.

OK. I see your point and fully agree with it!

[perazz]

@jvdp1 @jalvesz note that we already have many such constants (including precision-related) in the BLAS modules.

• For the user API, maybe all we need is a consistent way to export them? One way this was done in the codata PR was to use a mold parameter and functions
• For development, maybe the best approach would be to template them directly in fypp? this works well:

#! Return the unitary numeric constant for a given type and kind
#!
#! Args:
#!   type (string): an intrinsic type (complex, real, integer)
#!
#! Returns:
#!   Parameter value string
#!
#! E.g.,
#!   one('complex(dp)')  -> (1.0_dp, 0.0_dp)
#!   one('real(sp)') -> 1.0_sp
#!   one('integer') -> 1
#!
#:def one(type)
  #:assert 'integer' in type or 'real' in type or 'complex' in type
  #:if '(' in type and ')' in type
    #:set kind="_" + type.split('(')[1].split(')')[0]
  #:else
    #:set kind=""
  #:endif
  #:if 'integer' in type
    ${"1"+kind}$
  #:elif 'real' in type
    ${"1.0"+kind}$
  #:elif 'complex' in type
    ${"(1.0"+kind+",0.0"+kind+')'}$
  #:endif
#:enddef