linalg: determinant

[perazz]

Introduce determinant operator, based on LAPACK *GETRF functions.

• implementation
• tests
• exclude unsupported xdp
• documentation
• submodule

Prior art:

• Numpy: det(a)
• Scipy: det(a, overwrite_a=False, check_finite=True)
• IMSL: DET(a)

Proposed implementation: generic interface

1. d = det(a) -> pure function, no additional inputs, a is unchanged
2. d = det(a,overwrite_a=.false.,err=state) -> not pure, option to overwrite_a (do not allocate a temporary), return error state handler
3. .det.a Fortran operator interface, similar to 1), pure, but with the operator method

Error handling returns:

• LINALG_ERROR on singular matrix
• LINALG_VALUE_ERROR on invalid matrix size

cc @jvdp1 @jalvesz @fortran-lang/stdlib

[jvdp1]

Could it be a submdodule of stdblib_linalg, such that users can use it as:

use stdlib_linalg, only: det

instead of:

use stdlib_linalg_determinant, only: det

This will avoid an enormous amount of stdlib_linalg modules.
Inconvenient: this approach will require the compilation of the a large module.

What do you think? What would be the best strategy for the users?

[jalvesz]

This seems like a good approach, everything would fit together smoothly. This might require indeed a large header module but for the end user and even a developer who would like to pick just one method and work on it, it would be easier... I think.

[perazz]

Thanks @jvdp1 @jalvesz, great to iterate on the format, the more time we spend on this first LAPACK-backed function, the easier will be later.
The submodule one is a great idea and here is an implementation.

gfortran returns an error though:

Fixed thanks to Fortran Discourse

[jvdp1]

Thank you @perazz . Here are some suggestions. Overall, I like the API.

[jvdp1]

This let me think that we should review the structure of the specs. Currently we have 1 page per module. However, if these become available through stdlib_linalg, the specs should be modified accordingly (probably in another PR; I can start to work on a proposition when this is ready).

[jvdp1]

For future developers, it might be good to add a comment to explain why xdp is excluded.

[perazz]

504d90d

[perazz]

Even though it is not automated generation, should we manually implement blas and lapack backends for xdp? They would be the same as for quad-precision, because all constants are generalized, and 64-bit double precision is still the "lower precision" datatype for mixed-precision algorithms. We should just decide about real/complex prefixes, like s/c vs. d/z vs. q/w, for xdp they may be x/y?

[jvdp1]

Supporting xdp would be a nice feature IMO (even if I never use xdp). x/y prefixes make sense to me. If you think it is easy to implement it, I would support it. However, I don't think it is a priority right now.

[jalvesz]

If it serves for the discussion, I'm using the following:
in the common.fypp I added the notion of SUFFIX like

...
#! Collected (kind, type, suffix) tuples for real types
#:set REAL_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_KINDS))

#! Complex types to be considered during templating
#:set CMPLX_SUFFIX = ["c{}".format(k) for k in CMPLX_KINDS]

#! Collected (kind, type, suffix) tuples for complex types
#:set CMPLX_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))

Which then allows me to do this in the .fypp

#:for k1, t1, s1 in (KINDS_TYPES)
    type, public, extends(COO_t) :: COO_${s1}$
    ${t1}$, allocatable :: data(:) 
contains
     procedure :: get => get_value_coo_${s1}$
     procedure :: set => set_value_coo_${s1}$
end type
#:endfor

and obtaining this in the fortran file:

    type, public, extends(COO_t) :: COO_sp
        real(sp), allocatable :: data(:) 
    contains
        procedure :: get => get_value_coo_sp
        procedure :: set => set_value_coo_sp
    end type
    type, public, extends(COO_t) :: COO_dp
        real(dp), allocatable :: data(:) 
    contains
        procedure :: get => get_value_coo_dp
        procedure :: set => set_value_coo_dp
    end type
    type, public, extends(COO_t) :: COO_csp
        complex(sp), allocatable :: data(:) 
    contains
        procedure :: get => get_value_coo_csp
        procedure :: set => set_value_coo_csp
    end type
    type, public, extends(COO_t) :: COO_cdp
        complex(dp), allocatable :: data(:) 
    contains
        procedure :: get => get_value_coo_cdp
        procedure :: set => set_value_coo_cdp
    end type

So basically considering the real types as reference, and the complex types as extensions adding the c in front of the suffix.

[jvdp1]

Thank you @perazz . This PR is almost ready to be merged IMO. Specs should be still added.

[jvdp1]

Thank you @perazz. The specs LGTM. I proposed some minor suggestions.

[perazz]

Thanks a lot @jvdp1 - your edits look good! I've just added that it works for either real or complex matrices (it would seem it was for real only).

[perazz]

Thank you @jvdp1 @jalvesz for the reviews, let's wait a bit more, and if there are no further comments, I believe we could merge this before the end of the week.

[jvdp1]

I'll merge this PR. Thank you @perazz for this nice feature.